Loading...
Searching...
No Matches
mpr_inout.h File Reference

Go to the source code of this file.

## Macros

#define | DEFAULT_DIGITS 30
#define | MPR_DENSE 1
#define | MPR_SPARSE 2

## Functions

BOOLEAN | nuUResSolve (leftv res, leftv args)
| solve a multipolynomial system using the u-resultant Input ideal must be 0-dimensional and (currRing->N) == IDELEMS(ideal).
BOOLEAN | nuMPResMat (leftv res, leftv arg1, leftv arg2)
| returns module representing the multipolynomial resultant matrix Arguments 2: ideal i, int k k=0: use sparse resultant matrix of Gelfand, Kapranov and Zelevinsky k=1: use resultant matrix of Macaulay (k=0 is default)
BOOLEAN | nuLagSolve (leftv res, leftv arg1, leftv arg2, leftv arg3)
| find the (complex) roots an univariate polynomial Determines the roots of an univariate polynomial using Laguerres' root-solver.
BOOLEAN | nuVanderSys (leftv res, leftv arg1, leftv arg2, leftv arg3)
| COMPUTE: polynomial p with values given by v at points p1,..,pN derived from p; more precisely: consider p as point in K^n and v as N elements in K, let p1,..,pN be the points in K^n obtained by evaluating all monomials of degree 0,1,...,N at p in lexicographical order, then the procedure computes the polynomial f satisfying f(pi) = v[i] RETURN: polynomial f of degree d.
BOOLEAN | loNewtonP (leftv res, leftv arg1)
| compute Newton Polytopes of input polynomials
BOOLEAN | loSimplex (leftv res, leftv args)
| Implementation of the Simplex Algorithm.

## Macro Definition Documentation

## DEFAULT_DIGITS

#define DEFAULT_DIGITS 30

Definition at line 13 of file mpr_inout.h.

## MPR_DENSE

#define MPR_DENSE 1

Definition at line 15 of file mpr_inout.h.

## MPR_SPARSE

#define MPR_SPARSE 2

Definition at line 16 of file mpr_inout.h.

## Function Documentation

## loNewtonP()

BOOLEAN loNewtonP | ( | leftv | res,
|  | leftv | arg1 )

compute Newton Polytopes of input polynomials

Definition at line 4566 of file ipshell.cc.

{
res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
return FALSE;
}
FALSE
#define FALSE
Definition auxiliary.h:97
sleftv::Data
void * Data()
Definition subexpr.cc:1192
res
CanonicalForm res
Definition facAbsFact.cc:60
loNewtonPolytope
ideal loNewtonPolytope(const ideal id)
Definition mpr_base.cc:3191

## loSimplex()

BOOLEAN loSimplex | ( | leftv | res,
|  | leftv | args )

Implementation of the Simplex Algorithm.

For args, see class simplex.

Definition at line 4572 of file ipshell.cc.

{
if ( !(rField_is_long_R(currRing)) )
{
WerrorS("Ground field not implemented!");
return TRUE;
}
simplex * LP;
matrix m;
leftv v= args;
if ( v->Typ() != MATRIX_CMD ) // 1: matrix
return TRUE;
else
m= (matrix)(v->CopyD());
LP = new simplex(MATROWS(m),MATCOLS(m));
LP->mapFromMatrix(m);
v= v->next;
if ( v->Typ() != INT_CMD ) // 2: m = number of constraints
return TRUE;
else
LP->m= (int)(long)(v->Data());
v= v->next;
if ( v->Typ() != INT_CMD ) // 3: n = number of variables
return TRUE;
else
LP->n= (int)(long)(v->Data());
v= v->next;
if ( v->Typ() != INT_CMD ) // 4: m1 = number of <= constraints
return TRUE;
else
LP->m1= (int)(long)(v->Data());
v= v->next;
if ( v->Typ() != INT_CMD ) // 5: m2 = number of >= constraints
return TRUE;
else
LP->m2= (int)(long)(v->Data());
v= v->next;
if ( v->Typ() != INT_CMD ) // 6: m3 = number of == constraints
return TRUE;
else
LP->m3= (int)(long)(v->Data());
#ifdef mprDEBUG_PROT
Print("m (constraints) %d\n",LP->m);
Print("n (columns) %d\n",LP->n);
Print("m1 (<=) %d\n",LP->m1);
Print("m2 (>=) %d\n",LP->m2);
Print("m3 (==) %d\n",LP->m3);
#endif
LP->compute();
lists lres= (lists)omAlloc( sizeof(slists) );
lres->Init( 6 );
lres->m[0].rtyp= MATRIX_CMD; // output matrix
lres->m[0].data=(void*)LP->mapToMatrix(m);
lres->m[1].rtyp= INT_CMD; // found a solution?
lres->m[1].data=(void*)(long)LP->icase;
lres->m[2].rtyp= INTVEC_CMD;
lres->m[2].data=(void*)LP->posvToIV();
lres->m[3].rtyp= INTVEC_CMD;
lres->m[3].data=(void*)LP->zrovToIV();
lres->m[4].rtyp= INT_CMD;
lres->m[4].data=(void*)(long)LP->m;
lres->m[5].rtyp= INT_CMD;
lres->m[5].data=(void*)(long)LP->n;
res->data= (void*)lres;
return FALSE;
}
TRUE
#define TRUE
Definition auxiliary.h:101
m
int m
Definition cfEzgcd.cc:128
simplex
Linear Programming / Linear Optimization using Simplex - Algorithm.
Definition mpr_numeric.h:195
simplex::zrovToIV
intvec * zrovToIV()
Definition mpr_numeric.cc:1084
simplex::m1
int m1
Definition mpr_numeric.h:200
simplex::n
int n
Definition mpr_numeric.h:199
simplex::m2
int m2
Definition mpr_numeric.h:200
simplex::mapFromMatrix
BOOLEAN mapFromMatrix(matrix m)
Definition mpr_numeric.cc:1011
simplex::icase
int icase
Definition mpr_numeric.h:201
simplex::compute
void compute()
Definition mpr_numeric.cc:1095
simplex::m3
int m3
Definition mpr_numeric.h:200
simplex::m
int m
Definition mpr_numeric.h:198
simplex::mapToMatrix
matrix mapToMatrix(matrix m)
Definition mpr_numeric.cc:1040
simplex::posvToIV
intvec * posvToIV()
Definition mpr_numeric.cc:1073
sleftv::rtyp
int rtyp
Definition subexpr.h:91
sleftv::data
void * data
Definition subexpr.h:88
slists
Definition lists.h:24
slists::m
sleftv * m
Definition lists.h:46
slists::Init
INLINE_THIS void Init(int l=0)
Print
#define Print
Definition emacs.cc:80
v
const Variable & v
< [in] a sqrfree bivariate poly
Definition facBivar.h:39
WerrorS
void WerrorS(const char *s)
Definition feFopen.cc:24
MATRIX_CMD
@ MATRIX_CMD
Definition grammar.cc:287
matrix
ip_smatrix * matrix
Definition matpol.h:43
MATROWS
#define MATROWS(i)
Definition matpol.h:26
MATCOLS
#define MATCOLS(i)
Definition matpol.h:27
lists
slists * lists
Definition mpr_numeric.h:146
omAlloc
#define omAlloc(size)
Definition omAllocDecl.h:210
currRing
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition polys.cc:13
rField_is_long_R
static BOOLEAN rField_is_long_R(const ring r)
Definition ring.h:548
leftv
sleftv * leftv
Definition structs.h:53
INTVEC_CMD
@ INTVEC_CMD
Definition tok.h:101
INT_CMD
@ INT_CMD
Definition tok.h:96

## nuLagSolve()

BOOLEAN nuLagSolve | ( | leftv | res,
|  | leftv | arg1,
|  | leftv | arg2,
|  | leftv | arg3 )

find the (complex) roots an univariate polynomial Determines the roots of an univariate polynomial using Laguerres' root-solver.

Good for polynomials with low and middle degree (<40). Arguments 3: poly arg1 , int arg2 , int arg3 arg2>0: defines precision of fractional part if ground field is Q arg3: number of iterations for approximation of roots (default=2) Returns a list of all (complex) roots of the polynomial arg1

Definition at line 4681 of file ipshell.cc.

{
poly gls;
gls= (poly)(arg1->Data());
int howclean= (int)(long)arg3->Data();
if ( gls == NULL || pIsConstant( gls ) )
{
WerrorS("Input polynomial is constant!");
return TRUE;
}
if (rField_is_Zp(currRing))
{
int* r=Zp_roots(gls, currRing);
lists rlist;
rlist= (lists)omAlloc( sizeof(slists) );
rlist->Init( r[0] );
for(int i=r[0];i>0;i--)
{
rlist->m[i-1].data=n_Init(r[i],currRing->cf);
rlist->m[i-1].rtyp=NUMBER_CMD;
}
omFree(r);
res->data=rlist;
res->rtyp= LIST_CMD;
return FALSE;
}
if ( !(rField_is_R(currRing) ||
rField_is_Q(currRing) ||
rField_is_long_R(currRing) ||
rField_is_long_C(currRing)) )
{
WerrorS("Ground field not implemented!");
return TRUE;
}
if ( !(rField_is_R(currRing) || rField_is_long_R(currRing) || \
rField_is_long_C(currRing)) )
{
unsigned long int ii = (unsigned long int)arg2->Data();
setGMPFloatDigits( ii, ii );
}
int ldummy;
int deg= currRing->pLDeg( gls, &ldummy, currRing );
int i,vpos=0;
poly piter;
lists elist;
elist= (lists)omAlloc( sizeof(slists) );
elist->Init( 0 );
if ( rVar(currRing) > 1 )
{
piter= gls;
for ( i= 1; i <= rVar(currRing); i++ )
if ( pGetExp( piter, i ) )
{
vpos= i;
break;
}
while ( piter )
{
for ( i= 1; i <= rVar(currRing); i++ )
if ( (vpos != i) && (pGetExp( piter, i ) != 0) )
{
WerrorS("The input polynomial must be univariate!");
return TRUE;
}
pIter( piter );
}
}
rootContainer * roots= new rootContainer();
number * pcoeffs= (number *)omAlloc( (deg+1) * sizeof( number ) );
piter= gls;
for ( i= deg; i >= 0; i-- )
{
if ( piter && pTotaldegree(piter) == i )
{
pcoeffs[i]= nCopy( pGetCoeff( piter ) );
//nPrint( pcoeffs[i] );PrintS(" ");
pIter( piter );
}
else
{
pcoeffs[i]= nInit(0);
}
}
#ifdef mprDEBUG_PROT
for (i=deg; i >= 0; i--)
{
nPrint( pcoeffs[i] );PrintS(" ");
}
PrintLn();
#endif
roots->fillContainer( pcoeffs, NULL, 1, deg, rootContainer::onepoly, 1 );
roots->solver( howclean );
int elem= roots->getAnzRoots();
char *dummy;
int j;
lists rlist;
rlist= (lists)omAlloc( sizeof(slists) );
rlist->Init( elem );
if (rField_is_long_C(currRing))
{
for ( j= 0; j < elem; j++ )
{
rlist->m[j].rtyp=NUMBER_CMD;
rlist->m[j].data=(void *)nCopy((number)(roots->getRoot(j)));
//rlist->m[j].data=(void *)(number)(roots->getRoot(j));
}
}
else
{
for ( j= 0; j < elem; j++ )
{
dummy = complexToStr( (*roots)[j], gmp_output_digits, currRing->cf );
rlist->m[j].rtyp=STRING_CMD;
rlist->m[j].data=(void *)dummy;
}
}
elist->Clean();
//omFreeSize( (ADDRESS) elist, sizeof(slists) );
// this is (via fillContainer) the same data as in root
//for ( i= deg; i >= 0; i-- ) nDelete( &pcoeffs[i] );
//omFreeSize( (ADDRESS) pcoeffs, (deg+1) * sizeof( number ) );
delete roots;
res->data= (void*)rlist;
return FALSE;
}
i
int i
Definition cfEzgcd.cc:132
Zp_roots
int * Zp_roots(poly p, const ring r)
Definition clapsing.cc:2242
rootContainer
complex root finder for univariate polynomials based on laguers algorithm
Definition mpr_numeric.h:66
rootContainer::getRoot
gmp_complex * getRoot(const int i)
Definition mpr_numeric.h:88
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::onepoly
@ onepoly
Definition mpr_numeric.h:68
rootContainer::getAnzRoots
int getAnzRoots()
Definition mpr_numeric.h:97
rootContainer::solver
bool solver(const int polishmode=PM_NONE)
Definition mpr_numeric.cc:437
slists::Clean
void Clean(ring r=currRing)
Definition lists.h:26
n_Init
static FORCE_INLINE number n_Init(long i, const coeffs r)
a number representing i in the given coeff field/ring r
Definition coeffs.h:541
j
int j
Definition facHensel.cc:110
NUMBER_CMD
@ NUMBER_CMD
Definition grammar.cc:289
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
gmp_output_digits
EXTERN_VAR size_t gmp_output_digits
Definition mpr_base.h:115
complexToStr
char * complexToStr(gmp_complex &c, const unsigned int oprec, const coeffs src)
Definition mpr_complex.cc:693
setGMPFloatDigits
void setGMPFloatDigits(size_t digits, size_t rest)
Set size of mantissa digits - the number of output digits (basis 10) the size of mantissa consists of...
Definition mpr_complex.cc:60
nCopy
#define nCopy(n)
Definition numbers.h:15
nPrint
#define nPrint(a)
only for debug, over any initialized currRing
Definition numbers.h:46
nInit
#define nInit(i)
Definition numbers.h:24
omFree
#define omFree(addr)
Definition omAllocDecl.h:261
NULL
#define NULL
Definition omList.c:12
pTotaldegree
static long pTotaldegree(poly p)
Definition polys.h:283
pIsConstant
#define pIsConstant(p)
like above, except that Comp must be 0
Definition polys.h:239
pGetExp
#define pGetExp(p, i)
Exponent.
Definition polys.h:42
PrintS
void PrintS(const char *s)
Definition reporter.cc:288
PrintLn
void PrintLn()
Definition reporter.cc:314
rField_is_R
static BOOLEAN rField_is_R(const ring r)
Definition ring.h:524
rField_is_Zp
static BOOLEAN rField_is_Zp(const ring r)
Definition ring.h:506
rField_is_long_C
static BOOLEAN rField_is_long_C(const ring r)
Definition ring.h:551
rField_is_Q
static BOOLEAN rField_is_Q(const ring r)
Definition ring.h:512
rVar
static short rVar(const ring r)
define rVar(r) (r->N)
Definition ring.h:598
LIST_CMD
@ LIST_CMD
Definition tok.h:118
STRING_CMD
@ STRING_CMD
Definition tok.h:187

## nuMPResMat()

BOOLEAN nuMPResMat | ( | leftv | res,
|  | leftv | arg1,
|  | leftv | arg2 )

returns module representing the multipolynomial resultant matrix Arguments 2: ideal i, int k k=0: use sparse resultant matrix of Gelfand, Kapranov and Zelevinsky k=1: use resultant matrix of Macaulay (k=0 is default)

Definition at line 4658 of file ipshell.cc.

{
ideal gls = (ideal)(arg1->Data());
int imtype= (int)(long)arg2->Data();
uResultant::resMatType mtype= determineMType( imtype );
// check input ideal ( = polynomial system )
if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk )
{
return TRUE;
}
uResultant *resMat= new uResultant( gls, mtype, false );
if (resMat!=NULL)
{
res->rtyp = MODUL_CMD;
res->data= (void*)resMat->accessResMat()->getMatrix();
if (!errorreported) delete resMat;
}
return errorreported;
}
resMatrixBase::getMatrix
virtual ideal getMatrix()
Definition mpr_base.h:31
sleftv::Name
const char * Name()
Definition subexpr.h:120
uResultant
Base class for solving 0-dim poly systems using u-resultant.
Definition mpr_base.h:63
uResultant::resMatType
resMatType
Definition mpr_base.h:65
uResultant::accessResMat
resMatrixBase * accessResMat()
Definition mpr_base.h:78
errorreported
VAR short errorreported
Definition feFopen.cc:23
MODUL_CMD
@ MODUL_CMD
Definition grammar.cc:288
mprOk
@ mprOk
Definition mpr_base.h:98
determineMType
uResultant::resMatType determineMType(int imtype)
mprIdealCheck
mprState mprIdealCheck(const ideal theIdeal, const char *name, uResultant::resMatType mtype, BOOLEAN rmatrix=false)

## nuUResSolve()

BOOLEAN nuUResSolve | ( | leftv | res,
|  | leftv | args )

solve a multipolynomial system using the u-resultant Input ideal must be 0-dimensional and (currRing->N) == IDELEMS(ideal).

Resultant method can be MPR_DENSE, which uses Macaulay Resultant (good for dense homogeneous polynoms) or MPR_SPARSE, which uses Sparse Resultant (Gelfand, Kapranov, Zelevinsky). Arguments 4: ideal i, int k, int l, int m k=0: use sparse resultant matrix of Gelfand, Kapranov and Zelevinsky k=1: use resultant matrix of Macaulay (k=0 is default) l>0: defines precision of fractional part if ground field is Q m=0,1,2: number of iterations for approximation of roots (default=2) Returns a list containing the roots of the system.

Definition at line 4925 of file ipshell.cc.

{
leftv v= args;
ideal gls;
int imtype;
int howclean;
// get ideal
if ( v->Typ() != IDEAL_CMD )
return TRUE;
else gls= (ideal)(v->Data());
v= v->next;
// get resultant matrix type to use (0,1)
if ( v->Typ() != INT_CMD )
return TRUE;
else imtype= (int)(long)v->Data();
v= v->next;
if (imtype==0)
{
ideal test_id=idInit(1,1);
int j;
for(j=IDELEMS(gls)-1;j>=0;j--)
{
if (gls->m[j]!=NULL)
{
test_id->m[0]=gls->m[j];
intvec *dummy_w=id_QHomWeight(test_id, currRing);
if (dummy_w!=NULL)
{
WerrorS("Newton polytope not of expected dimension");
delete dummy_w;
return TRUE;
}
}
}
}
// get and set precision in digits ( > 0 )
if ( v->Typ() != INT_CMD )
return TRUE;
else if ( !(rField_is_R(currRing) || rField_is_long_R(currRing) || \
rField_is_long_C(currRing)) )
{
unsigned long int ii=(unsigned long int)v->Data();
setGMPFloatDigits( ii, ii );
}
v= v->next;
// get interpolation steps (0,1,2)
if ( v->Typ() != INT_CMD )
return TRUE;
else howclean= (int)(long)v->Data();
uResultant::resMatType mtype= determineMType( imtype );
int i,count;
lists listofroots= NULL;
number smv= NULL;
BOOLEAN interpolate_det= (mtype==uResultant::denseResMat)?TRUE:FALSE;
//emptylist= (lists)omAlloc( sizeof(slists) );
//emptylist->Init( 0 );
//res->rtyp = LIST_CMD;
//res->data= (void *)emptylist;
// check input ideal ( = polynomial system )
if ( mprIdealCheck( gls, args->Name(), mtype ) != mprOk )
{
return TRUE;
}
uResultant * ures;
rootContainer ** iproots;
rootContainer ** muiproots;
rootArranger * arranger;
// main task 1: setup of resultant matrix
ures= new uResultant( gls, mtype );
if ( ures->accessResMat()->initState() != resMatrixBase::ready )
{
WerrorS("Error occurred during matrix setup!");
return TRUE;
}
// if dense resultant, check if minor nonsingular
if ( mtype == uResultant::denseResMat )
{
smv= ures->accessResMat()->getSubDet();
#ifdef mprDEBUG_PROT
PrintS("// Determinant of submatrix: ");nPrint(smv);PrintLn();
#endif
if ( nIsZero(smv) )
{
WerrorS("Unsuitable input ideal: Minor of resultant matrix is singular!");
return TRUE;
}
}
// main task 2: Interpolate specialized resultant polynomials
if ( interpolate_det )
iproots= ures->interpolateDenseSP( false, smv );
else
iproots= ures->specializeInU( false, smv );
// main task 3: Interpolate specialized resultant polynomials
if ( interpolate_det )
muiproots= ures->interpolateDenseSP( true, smv );
else
muiproots= ures->specializeInU( true, smv );
#ifdef mprDEBUG_PROT
int c= iproots[0]->getAnzElems();
for (i=0; i < c; i++) pWrite(iproots[i]->getPoly());
c= muiproots[0]->getAnzElems();
for (i=0; i < c; i++) pWrite(muiproots[i]->getPoly());
#endif
// main task 4: Compute roots of specialized polys and match them up
arranger= new rootArranger( iproots, muiproots, howclean );
arranger->solve_all();
// get list of roots
if ( arranger->success() )
{
arranger->arrange();
listofroots= listOfRoots(arranger, gmp_output_digits );
}
else
{
WerrorS("Solver was unable to find any roots!");
return TRUE;
}
// free everything
count= iproots[0]->getAnzElems();
for (i=0; i < count; i++) delete iproots[i];
omFreeSize( (ADDRESS) iproots, count * sizeof(rootContainer*) );
count= muiproots[0]->getAnzElems();
for (i=0; i < count; i++) delete muiproots[i];
omFreeSize( (ADDRESS) muiproots, count * sizeof(rootContainer*) );
delete ures;
delete arranger;
if (smv!=NULL) nDelete( &smv );
res->data= (void *)listofroots;
//emptylist->Clean();
// omFreeSize( (ADDRESS) emptylist, sizeof(slists) );
return FALSE;
}
BOOLEAN
int BOOLEAN
Definition auxiliary.h:88
ADDRESS
void * ADDRESS
Definition auxiliary.h:120
intvec
Definition intvec.h:24
resMatrixBase::getSubDet
virtual number getSubDet()
Definition mpr_base.h:37
resMatrixBase::ready
@ ready
Definition mpr_base.h:26
resMatrixBase::initState
virtual IStateType initState() const
Definition mpr_base.h:41
rootArranger
Definition mpr_numeric.h:150
rootArranger::solve_all
void solve_all()
Definition mpr_numeric.cc:858
rootArranger::success
bool success()
Definition mpr_numeric.h:162
rootArranger::arrange
void arrange()
Definition mpr_numeric.cc:883
rootContainer::getAnzElems
int getAnzElems()
Definition mpr_numeric.h:95
uResultant::specializeInU
rootContainer ** specializeInU(BOOLEAN matchUp=false, const number subDetVal=NULL)
Definition mpr_base.cc:3060
uResultant::interpolateDenseSP
rootContainer ** interpolateDenseSP(BOOLEAN matchUp=false, const number subDetVal=NULL)
Definition mpr_base.cc:2922
uResultant::denseResMat
@ denseResMat
Definition mpr_base.h:65
IDEAL_CMD
@ IDEAL_CMD
Definition grammar.cc:285
listOfRoots
lists listOfRoots(rootArranger *self, const unsigned int oprec)
Definition ipshell.cc:5082
nDelete
#define nDelete(n)
Definition numbers.h:16
nIsZero
#define nIsZero(n)
Definition numbers.h:19
omFreeSize
#define omFreeSize(addr, size)
Definition omAllocDecl.h:260
pWrite
void pWrite(poly p)
Definition polys.h:309
count
int status int void size_t count
Definition si_signals.h:69
idInit
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition simpleideals.cc:35
id_QHomWeight
intvec * id_QHomWeight(ideal id, const ring r)
Definition simpleideals.cc:1917
IDELEMS
#define IDELEMS(i)
Definition simpleideals.h:23

## nuVanderSys()

BOOLEAN nuVanderSys | ( | leftv | res,
|  | leftv | arg1,
|  | leftv | arg2,
|  | leftv | arg3 )

COMPUTE: polynomial p with values given by v at points p1,..,pN derived from p; more precisely: consider p as point in K^n and v as N elements in K, let p1,..,pN be the points in K^n obtained by evaluating all monomials of degree 0,1,...,N at p in lexicographical order, then the procedure computes the polynomial f satisfying f(pi) = v[i] RETURN: polynomial f of degree d.

Definition at line 4824 of file ipshell.cc.

{
int i;
ideal p,w;
p= (ideal)arg1->Data();
w= (ideal)arg2->Data();
// w[0] = f(p^0)
// w[1] = f(p^1)
// ...
// p can be a vector of numbers (multivariate polynom)
// or one number (univariate polynom)
// tdg = deg(f)
int n= IDELEMS( p );
int m= IDELEMS( w );
int tdg= (int)(long)arg3->Data();
res->data= (void*)NULL;
// check the input
if ( tdg < 1 )
{
WerrorS("Last input parameter must be > 0!");
return TRUE;
}
if ( n != rVar(currRing) )
{
Werror("Size of first input ideal must be equal to %d!",rVar(currRing));
return TRUE;
}
if ( m != (int)pow((double)tdg+1,(double)n) )
{
Werror("Size of second input ideal must be equal to %d!",
(int)pow((double)tdg+1,(double)n));
return TRUE;
}
if ( !(rField_is_Q(currRing) /* ||
rField_is_R() || rField_is_long_R() ||
rField_is_long_C()*/ ) )
{
WerrorS("Ground field not implemented!");
return TRUE;
}
number tmp;
number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
for ( i= 0; i < n; i++ )
{
pevpoint[i]=nInit(0);
if ( (p->m)[i] )
{
tmp = pGetCoeff( (p->m)[i] );
if ( nIsZero(tmp) || nIsOne(tmp) || nIsMOne(tmp) )
{
omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
WerrorS("Elements of first input ideal must not be equal to -1, 0, 1!");
return TRUE;
}
} else tmp= NULL;
if ( !nIsZero(tmp) )
{
if ( !pIsConstant((p->m)[i]))
{
omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
WerrorS("Elements of first input ideal must be numbers!");
return TRUE;
}
pevpoint[i]= nCopy( tmp );
}
}
number *wresults= (number *)omAlloc( m * sizeof( number ) );
for ( i= 0; i < m; i++ )
{
wresults[i]= nInit(0);
if ( (w->m)[i] && !nIsZero(pGetCoeff((w->m)[i])) )
{
if ( !pIsConstant((w->m)[i]))
{
omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
WerrorS("Elements of second input ideal must be numbers!");
return TRUE;
}
wresults[i]= nCopy(pGetCoeff((w->m)[i]));
}
}
vandermonde vm( m, n, tdg, pevpoint, FALSE );
number *ncpoly= vm.interpolateDense( wresults );
// do not free ncpoly[]!!
poly rpoly= vm.numvec2poly( ncpoly );
omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
res->data= (void*)rpoly;
return FALSE;
}
pow
Rational pow(const Rational &a, int e)
Definition GMPrat.cc:411
p
int p
Definition cfModGcd.cc:4086
vandermonde
vandermonde system solver for interpolating polynomials from their values
Definition mpr_numeric.h:29
w
const CanonicalForm & w
Definition facAbsFact.cc:51
nIsMOne
#define nIsMOne(n)
Definition numbers.h:26
nIsOne
#define nIsOne(n)
Definition numbers.h:25
Werror
void Werror(const char *fmt,...)
Definition reporter.cc:189