My Project
Loading...
Searching...
No Matches
Macros | Functions
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). More...
 
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) More...
 
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. More...
 
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. More...
 
BOOLEAN loNewtonP (leftv res, leftv arg1)
 compute Newton Polytopes of input polynomials More...
 
BOOLEAN loSimplex (leftv res, leftv args)
 Implementation of the Simplex Algorithm. More...
 

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 4563 of file ipshell.cc.

4564{
4565 res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
4566 return FALSE;
4567}
#define FALSE
Definition: auxiliary.h:96
void * Data()
Definition: subexpr.cc:1162
CanonicalForm res
Definition: facAbsFact.cc:60
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 4569 of file ipshell.cc.

4570{
4571 if ( !(rField_is_long_R(currRing)) )
4572 {
4573 WerrorS("Ground field not implemented!");
4574 return TRUE;
4575 }
4576
4577 simplex * LP;
4578 matrix m;
4579
4580 leftv v= args;
4581 if ( v->Typ() != MATRIX_CMD ) // 1: matrix
4582 return TRUE;
4583 else
4584 m= (matrix)(v->CopyD());
4585
4586 LP = new simplex(MATROWS(m),MATCOLS(m));
4587 LP->mapFromMatrix(m);
4588
4589 v= v->next;
4590 if ( v->Typ() != INT_CMD ) // 2: m = number of constraints
4591 return TRUE;
4592 else
4593 LP->m= (int)(long)(v->Data());
4594
4595 v= v->next;
4596 if ( v->Typ() != INT_CMD ) // 3: n = number of variables
4597 return TRUE;
4598 else
4599 LP->n= (int)(long)(v->Data());
4600
4601 v= v->next;
4602 if ( v->Typ() != INT_CMD ) // 4: m1 = number of <= constraints
4603 return TRUE;
4604 else
4605 LP->m1= (int)(long)(v->Data());
4606
4607 v= v->next;
4608 if ( v->Typ() != INT_CMD ) // 5: m2 = number of >= constraints
4609 return TRUE;
4610 else
4611 LP->m2= (int)(long)(v->Data());
4612
4613 v= v->next;
4614 if ( v->Typ() != INT_CMD ) // 6: m3 = number of == constraints
4615 return TRUE;
4616 else
4617 LP->m3= (int)(long)(v->Data());
4618
4619#ifdef mprDEBUG_PROT
4620 Print("m (constraints) %d\n",LP->m);
4621 Print("n (columns) %d\n",LP->n);
4622 Print("m1 (<=) %d\n",LP->m1);
4623 Print("m2 (>=) %d\n",LP->m2);
4624 Print("m3 (==) %d\n",LP->m3);
4625#endif
4626
4627 LP->compute();
4628
4629 lists lres= (lists)omAlloc( sizeof(slists) );
4630 lres->Init( 6 );
4631
4632 lres->m[0].rtyp= MATRIX_CMD; // output matrix
4633 lres->m[0].data=(void*)LP->mapToMatrix(m);
4634
4635 lres->m[1].rtyp= INT_CMD; // found a solution?
4636 lres->m[1].data=(void*)(long)LP->icase;
4637
4638 lres->m[2].rtyp= INTVEC_CMD;
4639 lres->m[2].data=(void*)LP->posvToIV();
4640
4641 lres->m[3].rtyp= INTVEC_CMD;
4642 lres->m[3].data=(void*)LP->zrovToIV();
4643
4644 lres->m[4].rtyp= INT_CMD;
4645 lres->m[4].data=(void*)(long)LP->m;
4646
4647 lres->m[5].rtyp= INT_CMD;
4648 lres->m[5].data=(void*)(long)LP->n;
4649
4650 res->data= (void*)lres;
4651
4652 return FALSE;
4653}
#define TRUE
Definition: auxiliary.h:100
int m
Definition: cfEzgcd.cc:128
Variable next() const
Definition: variable.h:52
Linear Programming / Linear Optimization using Simplex - Algorithm.
Definition: mpr_numeric.h:195
intvec * zrovToIV()
BOOLEAN mapFromMatrix(matrix m)
int icase
Definition: mpr_numeric.h:201
void compute()
matrix mapToMatrix(matrix m)
intvec * posvToIV()
Class used for (list of) interpreter objects.
Definition: subexpr.h:83
int rtyp
Definition: subexpr.h:91
void * data
Definition: subexpr.h:88
Definition: lists.h:24
sleftv * m
Definition: lists.h:46
INLINE_THIS void Init(int l=0)
#define Print
Definition: emacs.cc:80
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
void WerrorS(const char *s)
Definition: feFopen.cc:24
@ MATRIX_CMD
Definition: grammar.cc:286
ip_smatrix * matrix
Definition: matpol.h:43
#define MATROWS(i)
Definition: matpol.h:26
#define MATCOLS(i)
Definition: matpol.h:27
slists * lists
Definition: mpr_numeric.h:146
#define omAlloc(size)
Definition: omAllocDecl.h:210
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
static BOOLEAN rField_is_long_R(const ring r)
Definition: ring.h:542
@ INTVEC_CMD
Definition: tok.h:101
@ 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 4678 of file ipshell.cc.

4679{
4680 poly gls;
4681 gls= (poly)(arg1->Data());
4682 int howclean= (int)(long)arg3->Data();
4683
4684 if ( gls == NULL || pIsConstant( gls ) )
4685 {
4686 WerrorS("Input polynomial is constant!");
4687 return TRUE;
4688 }
4689
4691 {
4692 int* r=Zp_roots(gls, currRing);
4693 lists rlist;
4694 rlist= (lists)omAlloc( sizeof(slists) );
4695 rlist->Init( r[0] );
4696 for(int i=r[0];i>0;i--)
4697 {
4698 rlist->m[i-1].data=n_Init(r[i],currRing->cf);
4699 rlist->m[i-1].rtyp=NUMBER_CMD;
4700 }
4701 omFree(r);
4702 res->data=rlist;
4703 res->rtyp= LIST_CMD;
4704 return FALSE;
4705 }
4706 if ( !(rField_is_R(currRing) ||
4710 {
4711 WerrorS("Ground field not implemented!");
4712 return TRUE;
4713 }
4714
4717 {
4718 unsigned long int ii = (unsigned long int)arg2->Data();
4719 setGMPFloatDigits( ii, ii );
4720 }
4721
4722 int ldummy;
4723 int deg= currRing->pLDeg( gls, &ldummy, currRing );
4724 int i,vpos=0;
4725 poly piter;
4726 lists elist;
4727
4728 elist= (lists)omAlloc( sizeof(slists) );
4729 elist->Init( 0 );
4730
4731 if ( rVar(currRing) > 1 )
4732 {
4733 piter= gls;
4734 for ( i= 1; i <= rVar(currRing); i++ )
4735 if ( pGetExp( piter, i ) )
4736 {
4737 vpos= i;
4738 break;
4739 }
4740 while ( piter )
4741 {
4742 for ( i= 1; i <= rVar(currRing); i++ )
4743 if ( (vpos != i) && (pGetExp( piter, i ) != 0) )
4744 {
4745 WerrorS("The input polynomial must be univariate!");
4746 return TRUE;
4747 }
4748 pIter( piter );
4749 }
4750 }
4751
4752 rootContainer * roots= new rootContainer();
4753 number * pcoeffs= (number *)omAlloc( (deg+1) * sizeof( number ) );
4754 piter= gls;
4755 for ( i= deg; i >= 0; i-- )
4756 {
4757 if ( piter && pTotaldegree(piter) == i )
4758 {
4759 pcoeffs[i]= nCopy( pGetCoeff( piter ) );
4760 //nPrint( pcoeffs[i] );PrintS(" ");
4761 pIter( piter );
4762 }
4763 else
4764 {
4765 pcoeffs[i]= nInit(0);
4766 }
4767 }
4768
4769#ifdef mprDEBUG_PROT
4770 for (i=deg; i >= 0; i--)
4771 {
4772 nPrint( pcoeffs[i] );PrintS(" ");
4773 }
4774 PrintLn();
4775#endif
4776
4777 roots->fillContainer( pcoeffs, NULL, 1, deg, rootContainer::onepoly, 1 );
4778 roots->solver( howclean );
4779
4780 int elem= roots->getAnzRoots();
4781 char *dummy;
4782 int j;
4783
4784 lists rlist;
4785 rlist= (lists)omAlloc( sizeof(slists) );
4786 rlist->Init( elem );
4787
4789 {
4790 for ( j= 0; j < elem; j++ )
4791 {
4792 rlist->m[j].rtyp=NUMBER_CMD;
4793 rlist->m[j].data=(void *)nCopy((number)(roots->getRoot(j)));
4794 //rlist->m[j].data=(void *)(number)(roots->getRoot(j));
4795 }
4796 }
4797 else
4798 {
4799 for ( j= 0; j < elem; j++ )
4800 {
4801 dummy = complexToStr( (*roots)[j], gmp_output_digits, currRing->cf );
4802 rlist->m[j].rtyp=STRING_CMD;
4803 rlist->m[j].data=(void *)dummy;
4804 }
4805 }
4806
4807 elist->Clean();
4808 //omFreeSize( (ADDRESS) elist, sizeof(slists) );
4809
4810 // this is (via fillContainer) the same data as in root
4811 //for ( i= deg; i >= 0; i-- ) nDelete( &pcoeffs[i] );
4812 //omFreeSize( (ADDRESS) pcoeffs, (deg+1) * sizeof( number ) );
4813
4814 delete roots;
4815
4816 res->data= (void*)rlist;
4817
4818 return FALSE;
4819}
int i
Definition: cfEzgcd.cc:132
int * Zp_roots(poly p, const ring r)
Definition: clapsing.cc:2188
complex root finder for univariate polynomials based on laguers algorithm
Definition: mpr_numeric.h:66
gmp_complex * getRoot(const int i)
Definition: mpr_numeric.h:88
void fillContainer(number *_coeffs, number *_ievpoint, const int _var, const int _tdg, const rootType _rt, const int _anz)
Definition: mpr_numeric.cc:300
int getAnzRoots()
Definition: mpr_numeric.h:97
bool solver(const int polishmode=PM_NONE)
Definition: mpr_numeric.cc:437
void Clean(ring r=currRing)
Definition: lists.h:26
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:535
int j
Definition: facHensel.cc:110
@ NUMBER_CMD
Definition: grammar.cc:288
#define pIter(p)
Definition: monomials.h:37
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
EXTERN_VAR size_t gmp_output_digits
Definition: mpr_base.h:115
char * complexToStr(gmp_complex &c, const unsigned int oprec, const coeffs src)
Definition: mpr_complex.cc:704
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
#define nCopy(n)
Definition: numbers.h:15
#define nPrint(a)
only for debug, over any initalized currRing
Definition: numbers.h:46
#define nInit(i)
Definition: numbers.h:24
#define omFree(addr)
Definition: omAllocDecl.h:261
#define NULL
Definition: omList.c:12
static long pTotaldegree(poly p)
Definition: polys.h:282
#define pIsConstant(p)
like above, except that Comp must be 0
Definition: polys.h:238
#define pGetExp(p, i)
Exponent.
Definition: polys.h:41
void PrintS(const char *s)
Definition: reporter.cc:284
void PrintLn()
Definition: reporter.cc:310
static BOOLEAN rField_is_R(const ring r)
Definition: ring.h:518
static BOOLEAN rField_is_Zp(const ring r)
Definition: ring.h:500
static BOOLEAN rField_is_long_C(const ring r)
Definition: ring.h:545
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:506
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:592
@ LIST_CMD
Definition: tok.h:118
@ STRING_CMD
Definition: tok.h:185

◆ 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 4655 of file ipshell.cc.

4656{
4657 ideal gls = (ideal)(arg1->Data());
4658 int imtype= (int)(long)arg2->Data();
4659
4660 uResultant::resMatType mtype= determineMType( imtype );
4661
4662 // check input ideal ( = polynomial system )
4663 if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk )
4664 {
4665 return TRUE;
4666 }
4667
4668 uResultant *resMat= new uResultant( gls, mtype, false );
4669 if (resMat!=NULL)
4670 {
4671 res->rtyp = MODUL_CMD;
4672 res->data= (void*)resMat->accessResMat()->getMatrix();
4673 if (!errorreported) delete resMat;
4674 }
4675 return errorreported;
4676}
virtual ideal getMatrix()
Definition: mpr_base.h:31
const char * Name()
Definition: subexpr.h:120
Base class for solving 0-dim poly systems using u-resultant.
Definition: mpr_base.h:63
resMatrixBase * accessResMat()
Definition: mpr_base.h:78
VAR short errorreported
Definition: feFopen.cc:23
@ MODUL_CMD
Definition: grammar.cc:287
@ mprOk
Definition: mpr_base.h:98
uResultant::resMatType determineMType(int imtype)
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 4922 of file ipshell.cc.

4923{
4924 leftv v= args;
4925
4926 ideal gls;
4927 int imtype;
4928 int howclean;
4929
4930 // get ideal
4931 if ( v->Typ() != IDEAL_CMD )
4932 return TRUE;
4933 else gls= (ideal)(v->Data());
4934 v= v->next;
4935
4936 // get resultant matrix type to use (0,1)
4937 if ( v->Typ() != INT_CMD )
4938 return TRUE;
4939 else imtype= (int)(long)v->Data();
4940 v= v->next;
4941
4942 if (imtype==0)
4943 {
4944 ideal test_id=idInit(1,1);
4945 int j;
4946 for(j=IDELEMS(gls)-1;j>=0;j--)
4947 {
4948 if (gls->m[j]!=NULL)
4949 {
4950 test_id->m[0]=gls->m[j];
4951 intvec *dummy_w=id_QHomWeight(test_id, currRing);
4952 if (dummy_w!=NULL)
4953 {
4954 WerrorS("Newton polytope not of expected dimension");
4955 delete dummy_w;
4956 return TRUE;
4957 }
4958 }
4959 }
4960 }
4961
4962 // get and set precision in digits ( > 0 )
4963 if ( v->Typ() != INT_CMD )
4964 return TRUE;
4965 else if ( !(rField_is_R(currRing) || rField_is_long_R(currRing) || \
4967 {
4968 unsigned long int ii=(unsigned long int)v->Data();
4969 setGMPFloatDigits( ii, ii );
4970 }
4971 v= v->next;
4972
4973 // get interpolation steps (0,1,2)
4974 if ( v->Typ() != INT_CMD )
4975 return TRUE;
4976 else howclean= (int)(long)v->Data();
4977
4978 uResultant::resMatType mtype= determineMType( imtype );
4979 int i,count;
4980 lists listofroots= NULL;
4981 number smv= NULL;
4982 BOOLEAN interpolate_det= (mtype==uResultant::denseResMat)?TRUE:FALSE;
4983
4984 //emptylist= (lists)omAlloc( sizeof(slists) );
4985 //emptylist->Init( 0 );
4986
4987 //res->rtyp = LIST_CMD;
4988 //res->data= (void *)emptylist;
4989
4990 // check input ideal ( = polynomial system )
4991 if ( mprIdealCheck( gls, args->Name(), mtype ) != mprOk )
4992 {
4993 return TRUE;
4994 }
4995
4996 uResultant * ures;
4997 rootContainer ** iproots;
4998 rootContainer ** muiproots;
4999 rootArranger * arranger;
5000
5001 // main task 1: setup of resultant matrix
5002 ures= new uResultant( gls, mtype );
5003 if ( ures->accessResMat()->initState() != resMatrixBase::ready )
5004 {
5005 WerrorS("Error occurred during matrix setup!");
5006 return TRUE;
5007 }
5008
5009 // if dense resultant, check if minor nonsingular
5010 if ( mtype == uResultant::denseResMat )
5011 {
5012 smv= ures->accessResMat()->getSubDet();
5013#ifdef mprDEBUG_PROT
5014 PrintS("// Determinant of submatrix: ");nPrint(smv);PrintLn();
5015#endif
5016 if ( nIsZero(smv) )
5017 {
5018 WerrorS("Unsuitable input ideal: Minor of resultant matrix is singular!");
5019 return TRUE;
5020 }
5021 }
5022
5023 // main task 2: Interpolate specialized resultant polynomials
5024 if ( interpolate_det )
5025 iproots= ures->interpolateDenseSP( false, smv );
5026 else
5027 iproots= ures->specializeInU( false, smv );
5028
5029 // main task 3: Interpolate specialized resultant polynomials
5030 if ( interpolate_det )
5031 muiproots= ures->interpolateDenseSP( true, smv );
5032 else
5033 muiproots= ures->specializeInU( true, smv );
5034
5035#ifdef mprDEBUG_PROT
5036 int c= iproots[0]->getAnzElems();
5037 for (i=0; i < c; i++) pWrite(iproots[i]->getPoly());
5038 c= muiproots[0]->getAnzElems();
5039 for (i=0; i < c; i++) pWrite(muiproots[i]->getPoly());
5040#endif
5041
5042 // main task 4: Compute roots of specialized polys and match them up
5043 arranger= new rootArranger( iproots, muiproots, howclean );
5044 arranger->solve_all();
5045
5046 // get list of roots
5047 if ( arranger->success() )
5048 {
5049 arranger->arrange();
5050 listofroots= listOfRoots(arranger, gmp_output_digits );
5051 }
5052 else
5053 {
5054 WerrorS("Solver was unable to find any roots!");
5055 return TRUE;
5056 }
5057
5058 // free everything
5059 count= iproots[0]->getAnzElems();
5060 for (i=0; i < count; i++) delete iproots[i];
5061 omFreeSize( (ADDRESS) iproots, count * sizeof(rootContainer*) );
5062 count= muiproots[0]->getAnzElems();
5063 for (i=0; i < count; i++) delete muiproots[i];
5064 omFreeSize( (ADDRESS) muiproots, count * sizeof(rootContainer*) );
5065
5066 delete ures;
5067 delete arranger;
5068 if (smv!=NULL) nDelete( &smv );
5069
5070 res->data= (void *)listofroots;
5071
5072 //emptylist->Clean();
5073 // omFreeSize( (ADDRESS) emptylist, sizeof(slists) );
5074
5075 return FALSE;
5076}
int BOOLEAN
Definition: auxiliary.h:87
void * ADDRESS
Definition: auxiliary.h:119
Definition: intvec.h:23
virtual number getSubDet()
Definition: mpr_base.h:37
virtual IStateType initState() const
Definition: mpr_base.h:41
void solve_all()
Definition: mpr_numeric.cc:858
bool success()
Definition: mpr_numeric.h:162
void arrange()
Definition: mpr_numeric.cc:883
int getAnzElems()
Definition: mpr_numeric.h:95
rootContainer ** specializeInU(BOOLEAN matchUp=false, const number subDetVal=NULL)
Definition: mpr_base.cc:3060
rootContainer ** interpolateDenseSP(BOOLEAN matchUp=false, const number subDetVal=NULL)
Definition: mpr_base.cc:2922
@ denseResMat
Definition: mpr_base.h:65
@ IDEAL_CMD
Definition: grammar.cc:284
lists listOfRoots(rootArranger *self, const unsigned int oprec)
Definition: ipshell.cc:5079
#define nDelete(n)
Definition: numbers.h:16
#define nIsZero(n)
Definition: numbers.h:19
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
void pWrite(poly p)
Definition: polys.h:308
int status int void size_t count
Definition: si_signals.h:59
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:35
intvec * id_QHomWeight(ideal id, const ring r)
#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 4821 of file ipshell.cc.

4822{
4823 int i;
4824 ideal p,w;
4825 p= (ideal)arg1->Data();
4826 w= (ideal)arg2->Data();
4827
4828 // w[0] = f(p^0)
4829 // w[1] = f(p^1)
4830 // ...
4831 // p can be a vector of numbers (multivariate polynom)
4832 // or one number (univariate polynom)
4833 // tdg = deg(f)
4834
4835 int n= IDELEMS( p );
4836 int m= IDELEMS( w );
4837 int tdg= (int)(long)arg3->Data();
4838
4839 res->data= (void*)NULL;
4840
4841 // check the input
4842 if ( tdg < 1 )
4843 {
4844 WerrorS("Last input parameter must be > 0!");
4845 return TRUE;
4846 }
4847 if ( n != rVar(currRing) )
4848 {
4849 Werror("Size of first input ideal must be equal to %d!",rVar(currRing));
4850 return TRUE;
4851 }
4852 if ( m != (int)pow((double)tdg+1,(double)n) )
4853 {
4854 Werror("Size of second input ideal must be equal to %d!",
4855 (int)pow((double)tdg+1,(double)n));
4856 return TRUE;
4857 }
4858 if ( !(rField_is_Q(currRing) /* ||
4859 rField_is_R() || rField_is_long_R() ||
4860 rField_is_long_C()*/ ) )
4861 {
4862 WerrorS("Ground field not implemented!");
4863 return TRUE;
4864 }
4865
4866 number tmp;
4867 number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
4868 for ( i= 0; i < n; i++ )
4869 {
4870 pevpoint[i]=nInit(0);
4871 if ( (p->m)[i] )
4872 {
4873 tmp = pGetCoeff( (p->m)[i] );
4874 if ( nIsZero(tmp) || nIsOne(tmp) || nIsMOne(tmp) )
4875 {
4876 omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4877 WerrorS("Elements of first input ideal must not be equal to -1, 0, 1!");
4878 return TRUE;
4879 }
4880 } else tmp= NULL;
4881 if ( !nIsZero(tmp) )
4882 {
4883 if ( !pIsConstant((p->m)[i]))
4884 {
4885 omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4886 WerrorS("Elements of first input ideal must be numbers!");
4887 return TRUE;
4888 }
4889 pevpoint[i]= nCopy( tmp );
4890 }
4891 }
4892
4893 number *wresults= (number *)omAlloc( m * sizeof( number ) );
4894 for ( i= 0; i < m; i++ )
4895 {
4896 wresults[i]= nInit(0);
4897 if ( (w->m)[i] && !nIsZero(pGetCoeff((w->m)[i])) )
4898 {
4899 if ( !pIsConstant((w->m)[i]))
4900 {
4901 omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4902 omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4903 WerrorS("Elements of second input ideal must be numbers!");
4904 return TRUE;
4905 }
4906 wresults[i]= nCopy(pGetCoeff((w->m)[i]));
4907 }
4908 }
4909
4910 vandermonde vm( m, n, tdg, pevpoint, FALSE );
4911 number *ncpoly= vm.interpolateDense( wresults );
4912 // do not free ncpoly[]!!
4913 poly rpoly= vm.numvec2poly( ncpoly );
4914
4915 omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4916 omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4917
4918 res->data= (void*)rpoly;
4919 return FALSE;
4920}
Rational pow(const Rational &a, int e)
Definition: GMPrat.cc:411
int p
Definition: cfModGcd.cc:4078
vandermonde system solver for interpolating polynomials from their values
Definition: mpr_numeric.h:29
const CanonicalForm & w
Definition: facAbsFact.cc:51
#define nIsMOne(n)
Definition: numbers.h:26
#define nIsOne(n)
Definition: numbers.h:25
void Werror(const char *fmt,...)
Definition: reporter.cc:189