My Project
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 4647 of file ipshell.cc.

4648 {
4649  res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
4650  return FALSE;
4651 }
#define FALSE
Definition: auxiliary.h:96
void * Data()
Definition: subexpr.cc:1154
CanonicalForm res
Definition: facAbsFact.cc:60
ideal loNewtonPolytope(const ideal id)
Definition: mpr_base.cc:3190

◆ loSimplex()

BOOLEAN loSimplex ( leftv  res,
leftv  args 
)

Implementation of the Simplex Algorithm.

For args, see class simplex.

Definition at line 4653 of file ipshell.cc.

4654 {
4655  if ( !(rField_is_long_R(currRing)) )
4656  {
4657  WerrorS("Ground field not implemented!");
4658  return TRUE;
4659  }
4660 
4661  simplex * LP;
4662  matrix m;
4663 
4664  leftv v= args;
4665  if ( v->Typ() != MATRIX_CMD ) // 1: matrix
4666  return TRUE;
4667  else
4668  m= (matrix)(v->CopyD());
4669 
4670  LP = new simplex(MATROWS(m),MATCOLS(m));
4671  LP->mapFromMatrix(m);
4672 
4673  v= v->next;
4674  if ( v->Typ() != INT_CMD ) // 2: m = number of constraints
4675  return TRUE;
4676  else
4677  LP->m= (int)(long)(v->Data());
4678 
4679  v= v->next;
4680  if ( v->Typ() != INT_CMD ) // 3: n = number of variables
4681  return TRUE;
4682  else
4683  LP->n= (int)(long)(v->Data());
4684 
4685  v= v->next;
4686  if ( v->Typ() != INT_CMD ) // 4: m1 = number of <= constraints
4687  return TRUE;
4688  else
4689  LP->m1= (int)(long)(v->Data());
4690 
4691  v= v->next;
4692  if ( v->Typ() != INT_CMD ) // 5: m2 = number of >= constraints
4693  return TRUE;
4694  else
4695  LP->m2= (int)(long)(v->Data());
4696 
4697  v= v->next;
4698  if ( v->Typ() != INT_CMD ) // 6: m3 = number of == constraints
4699  return TRUE;
4700  else
4701  LP->m3= (int)(long)(v->Data());
4702 
4703 #ifdef mprDEBUG_PROT
4704  Print("m (constraints) %d\n",LP->m);
4705  Print("n (columns) %d\n",LP->n);
4706  Print("m1 (<=) %d\n",LP->m1);
4707  Print("m2 (>=) %d\n",LP->m2);
4708  Print("m3 (==) %d\n",LP->m3);
4709 #endif
4710 
4711  LP->compute();
4712 
4713  lists lres= (lists)omAlloc( sizeof(slists) );
4714  lres->Init( 6 );
4715 
4716  lres->m[0].rtyp= MATRIX_CMD; // output matrix
4717  lres->m[0].data=(void*)LP->mapToMatrix(m);
4718 
4719  lres->m[1].rtyp= INT_CMD; // found a solution?
4720  lres->m[1].data=(void*)(long)LP->icase;
4721 
4722  lres->m[2].rtyp= INTVEC_CMD;
4723  lres->m[2].data=(void*)LP->posvToIV();
4724 
4725  lres->m[3].rtyp= INTVEC_CMD;
4726  lres->m[3].data=(void*)LP->zrovToIV();
4727 
4728  lres->m[4].rtyp= INT_CMD;
4729  lres->m[4].data=(void*)(long)LP->m;
4730 
4731  lres->m[5].rtyp= INT_CMD;
4732  lres->m[5].data=(void*)(long)LP->n;
4733 
4734  res->data= (void*)lres;
4735 
4736  return FALSE;
4737 }
#define TRUE
Definition: auxiliary.h:100
int m
Definition: cfEzgcd.cc:128
Variable next() const
Definition: factory.h:153
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:544
@ 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 4762 of file ipshell.cc.

4763 {
4764  poly gls;
4765  gls= (poly)(arg1->Data());
4766  int howclean= (int)(long)arg3->Data();
4767 
4768  if ( gls == NULL || pIsConstant( gls ) )
4769  {
4770  WerrorS("Input polynomial is constant!");
4771  return TRUE;
4772  }
4773 
4774  if (rField_is_Zp(currRing))
4775  {
4776  int* r=Zp_roots(gls, currRing);
4777  lists rlist;
4778  rlist= (lists)omAlloc( sizeof(slists) );
4779  rlist->Init( r[0] );
4780  for(int i=r[0];i>0;i--)
4781  {
4782  rlist->m[i-1].data=n_Init(r[i],currRing);
4783  rlist->m[i-1].rtyp=NUMBER_CMD;
4784  }
4785  omFree(r);
4786  res->data=rlist;
4787  res->rtyp= LIST_CMD;
4788  return FALSE;
4789  }
4790  if ( !(rField_is_R(currRing) ||
4791  rField_is_Q(currRing) ||
4794  {
4795  WerrorS("Ground field not implemented!");
4796  return TRUE;
4797  }
4798 
4801  {
4802  unsigned long int ii = (unsigned long int)arg2->Data();
4803  setGMPFloatDigits( ii, ii );
4804  }
4805 
4806  int ldummy;
4807  int deg= currRing->pLDeg( gls, &ldummy, currRing );
4808  int i,vpos=0;
4809  poly piter;
4810  lists elist;
4811 
4812  elist= (lists)omAlloc( sizeof(slists) );
4813  elist->Init( 0 );
4814 
4815  if ( rVar(currRing) > 1 )
4816  {
4817  piter= gls;
4818  for ( i= 1; i <= rVar(currRing); i++ )
4819  if ( pGetExp( piter, i ) )
4820  {
4821  vpos= i;
4822  break;
4823  }
4824  while ( piter )
4825  {
4826  for ( i= 1; i <= rVar(currRing); i++ )
4827  if ( (vpos != i) && (pGetExp( piter, i ) != 0) )
4828  {
4829  WerrorS("The input polynomial must be univariate!");
4830  return TRUE;
4831  }
4832  pIter( piter );
4833  }
4834  }
4835 
4836  rootContainer * roots= new rootContainer();
4837  number * pcoeffs= (number *)omAlloc( (deg+1) * sizeof( number ) );
4838  piter= gls;
4839  for ( i= deg; i >= 0; i-- )
4840  {
4841  if ( piter && pTotaldegree(piter) == i )
4842  {
4843  pcoeffs[i]= nCopy( pGetCoeff( piter ) );
4844  //nPrint( pcoeffs[i] );PrintS(" ");
4845  pIter( piter );
4846  }
4847  else
4848  {
4849  pcoeffs[i]= nInit(0);
4850  }
4851  }
4852 
4853 #ifdef mprDEBUG_PROT
4854  for (i=deg; i >= 0; i--)
4855  {
4856  nPrint( pcoeffs[i] );PrintS(" ");
4857  }
4858  PrintLn();
4859 #endif
4860 
4861  roots->fillContainer( pcoeffs, NULL, 1, deg, rootContainer::onepoly, 1 );
4862  roots->solver( howclean );
4863 
4864  int elem= roots->getAnzRoots();
4865  char *dummy;
4866  int j;
4867 
4868  lists rlist;
4869  rlist= (lists)omAlloc( sizeof(slists) );
4870  rlist->Init( elem );
4871 
4873  {
4874  for ( j= 0; j < elem; j++ )
4875  {
4876  rlist->m[j].rtyp=NUMBER_CMD;
4877  rlist->m[j].data=(void *)nCopy((number)(roots->getRoot(j)));
4878  //rlist->m[j].data=(void *)(number)(roots->getRoot(j));
4879  }
4880  }
4881  else
4882  {
4883  for ( j= 0; j < elem; j++ )
4884  {
4885  dummy = complexToStr( (*roots)[j], gmp_output_digits, currRing->cf );
4886  rlist->m[j].rtyp=STRING_CMD;
4887  rlist->m[j].data=(void *)dummy;
4888  }
4889  }
4890 
4891  elist->Clean();
4892  //omFreeSize( (ADDRESS) elist, sizeof(slists) );
4893 
4894  // this is (via fillContainer) the same data as in root
4895  //for ( i= deg; i >= 0; i-- ) nDelete( &pcoeffs[i] );
4896  //omFreeSize( (ADDRESS) pcoeffs, (deg+1) * sizeof( number ) );
4897 
4898  delete roots;
4899 
4900  res->data= (void*)rlist;
4901 
4902  return FALSE;
4903 }
int i
Definition: cfEzgcd.cc:132
int * Zp_roots(poly p, const ring r)
Definition: clapsing.cc:2048
complex root finder for univariate polynomials based on laguers algorithm
Definition: mpr_numeric.h:66
void fillContainer(number *_coeffs, number *_ievpoint, const int _var, const int _tdg, const rootType _rt, const int _anz)
Definition: mpr_numeric.cc:299
gmp_complex * getRoot(const int i)
Definition: mpr_numeric.h:88
int getAnzRoots()
Definition: mpr_numeric.h:97
bool solver(const int polishmode=PM_NONE)
Definition: mpr_numeric.cc:436
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:539
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:520
static BOOLEAN rField_is_Zp(const ring r)
Definition: ring.h:502
static BOOLEAN rField_is_long_C(const ring r)
Definition: ring.h:547
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:508
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:594
@ 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 4739 of file ipshell.cc.

4740 {
4741  ideal gls = (ideal)(arg1->Data());
4742  int imtype= (int)(long)arg2->Data();
4743 
4744  uResultant::resMatType mtype= determineMType( imtype );
4745 
4746  // check input ideal ( = polynomial system )
4747  if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk )
4748  {
4749  return TRUE;
4750  }
4751 
4752  uResultant *resMat= new uResultant( gls, mtype, false );
4753  if (resMat!=NULL)
4754  {
4755  res->rtyp = MODUL_CMD;
4756  res->data= (void*)resMat->accessResMat()->getMatrix();
4757  if (!errorreported) delete resMat;
4758  }
4759  return errorreported;
4760 }
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 5006 of file ipshell.cc.

5007 {
5008  leftv v= args;
5009 
5010  ideal gls;
5011  int imtype;
5012  int howclean;
5013 
5014  // get ideal
5015  if ( v->Typ() != IDEAL_CMD )
5016  return TRUE;
5017  else gls= (ideal)(v->Data());
5018  v= v->next;
5019 
5020  // get resultant matrix type to use (0,1)
5021  if ( v->Typ() != INT_CMD )
5022  return TRUE;
5023  else imtype= (int)(long)v->Data();
5024  v= v->next;
5025 
5026  if (imtype==0)
5027  {
5028  ideal test_id=idInit(1,1);
5029  int j;
5030  for(j=IDELEMS(gls)-1;j>=0;j--)
5031  {
5032  if (gls->m[j]!=NULL)
5033  {
5034  test_id->m[0]=gls->m[j];
5035  intvec *dummy_w=id_QHomWeight(test_id, currRing);
5036  if (dummy_w!=NULL)
5037  {
5038  WerrorS("Newton polytope not of expected dimension");
5039  delete dummy_w;
5040  return TRUE;
5041  }
5042  }
5043  }
5044  }
5045 
5046  // get and set precision in digits ( > 0 )
5047  if ( v->Typ() != INT_CMD )
5048  return TRUE;
5049  else if ( !(rField_is_R(currRing) || rField_is_long_R(currRing) || \
5051  {
5052  unsigned long int ii=(unsigned long int)v->Data();
5053  setGMPFloatDigits( ii, ii );
5054  }
5055  v= v->next;
5056 
5057  // get interpolation steps (0,1,2)
5058  if ( v->Typ() != INT_CMD )
5059  return TRUE;
5060  else howclean= (int)(long)v->Data();
5061 
5062  uResultant::resMatType mtype= determineMType( imtype );
5063  int i,count;
5064  lists listofroots= NULL;
5065  number smv= NULL;
5066  BOOLEAN interpolate_det= (mtype==uResultant::denseResMat)?TRUE:FALSE;
5067 
5068  //emptylist= (lists)omAlloc( sizeof(slists) );
5069  //emptylist->Init( 0 );
5070 
5071  //res->rtyp = LIST_CMD;
5072  //res->data= (void *)emptylist;
5073 
5074  // check input ideal ( = polynomial system )
5075  if ( mprIdealCheck( gls, args->Name(), mtype ) != mprOk )
5076  {
5077  return TRUE;
5078  }
5079 
5080  uResultant * ures;
5081  rootContainer ** iproots;
5082  rootContainer ** muiproots;
5083  rootArranger * arranger;
5084 
5085  // main task 1: setup of resultant matrix
5086  ures= new uResultant( gls, mtype );
5087  if ( ures->accessResMat()->initState() != resMatrixBase::ready )
5088  {
5089  WerrorS("Error occurred during matrix setup!");
5090  return TRUE;
5091  }
5092 
5093  // if dense resultant, check if minor nonsingular
5094  if ( mtype == uResultant::denseResMat )
5095  {
5096  smv= ures->accessResMat()->getSubDet();
5097 #ifdef mprDEBUG_PROT
5098  PrintS("// Determinant of submatrix: ");nPrint(smv);PrintLn();
5099 #endif
5100  if ( nIsZero(smv) )
5101  {
5102  WerrorS("Unsuitable input ideal: Minor of resultant matrix is singular!");
5103  return TRUE;
5104  }
5105  }
5106 
5107  // main task 2: Interpolate specialized resultant polynomials
5108  if ( interpolate_det )
5109  iproots= ures->interpolateDenseSP( false, smv );
5110  else
5111  iproots= ures->specializeInU( false, smv );
5112 
5113  // main task 3: Interpolate specialized resultant polynomials
5114  if ( interpolate_det )
5115  muiproots= ures->interpolateDenseSP( true, smv );
5116  else
5117  muiproots= ures->specializeInU( true, smv );
5118 
5119 #ifdef mprDEBUG_PROT
5120  int c= iproots[0]->getAnzElems();
5121  for (i=0; i < c; i++) pWrite(iproots[i]->getPoly());
5122  c= muiproots[0]->getAnzElems();
5123  for (i=0; i < c; i++) pWrite(muiproots[i]->getPoly());
5124 #endif
5125 
5126  // main task 4: Compute roots of specialized polys and match them up
5127  arranger= new rootArranger( iproots, muiproots, howclean );
5128  arranger->solve_all();
5129 
5130  // get list of roots
5131  if ( arranger->success() )
5132  {
5133  arranger->arrange();
5134  listofroots= listOfRoots(arranger, gmp_output_digits );
5135  }
5136  else
5137  {
5138  WerrorS("Solver was unable to find any roots!");
5139  return TRUE;
5140  }
5141 
5142  // free everything
5143  count= iproots[0]->getAnzElems();
5144  for (i=0; i < count; i++) delete iproots[i];
5145  omFreeSize( (ADDRESS) iproots, count * sizeof(rootContainer*) );
5146  count= muiproots[0]->getAnzElems();
5147  for (i=0; i < count; i++) delete muiproots[i];
5148  omFreeSize( (ADDRESS) muiproots, count * sizeof(rootContainer*) );
5149 
5150  delete ures;
5151  delete arranger;
5152  nDelete( &smv );
5153 
5154  res->data= (void *)listofroots;
5155 
5156  //emptylist->Clean();
5157  // omFreeSize( (ADDRESS) emptylist, sizeof(slists) );
5158 
5159  return FALSE;
5160 }
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:857
bool success()
Definition: mpr_numeric.h:162
void arrange()
Definition: mpr_numeric.cc:882
int getAnzElems()
Definition: mpr_numeric.h:95
rootContainer ** specializeInU(BOOLEAN matchUp=false, const number subDetVal=NULL)
Definition: mpr_base.cc:3059
rootContainer ** interpolateDenseSP(BOOLEAN matchUp=false, const number subDetVal=NULL)
Definition: mpr_base.cc:2921
@ denseResMat
Definition: mpr_base.h:65
@ IDEAL_CMD
Definition: grammar.cc:284
lists listOfRoots(rootArranger *self, const unsigned int oprec)
Definition: ipshell.cc:5163
#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 4905 of file ipshell.cc.

4906 {
4907  int i;
4908  ideal p,w;
4909  p= (ideal)arg1->Data();
4910  w= (ideal)arg2->Data();
4911 
4912  // w[0] = f(p^0)
4913  // w[1] = f(p^1)
4914  // ...
4915  // p can be a vector of numbers (multivariate polynom)
4916  // or one number (univariate polynom)
4917  // tdg = deg(f)
4918 
4919  int n= IDELEMS( p );
4920  int m= IDELEMS( w );
4921  int tdg= (int)(long)arg3->Data();
4922 
4923  res->data= (void*)NULL;
4924 
4925  // check the input
4926  if ( tdg < 1 )
4927  {
4928  WerrorS("Last input parameter must be > 0!");
4929  return TRUE;
4930  }
4931  if ( n != rVar(currRing) )
4932  {
4933  Werror("Size of first input ideal must be equal to %d!",rVar(currRing));
4934  return TRUE;
4935  }
4936  if ( m != (int)pow((double)tdg+1,(double)n) )
4937  {
4938  Werror("Size of second input ideal must be equal to %d!",
4939  (int)pow((double)tdg+1,(double)n));
4940  return TRUE;
4941  }
4942  if ( !(rField_is_Q(currRing) /* ||
4943  rField_is_R() || rField_is_long_R() ||
4944  rField_is_long_C()*/ ) )
4945  {
4946  WerrorS("Ground field not implemented!");
4947  return TRUE;
4948  }
4949 
4950  number tmp;
4951  number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
4952  for ( i= 0; i < n; i++ )
4953  {
4954  pevpoint[i]=nInit(0);
4955  if ( (p->m)[i] )
4956  {
4957  tmp = pGetCoeff( (p->m)[i] );
4958  if ( nIsZero(tmp) || nIsOne(tmp) || nIsMOne(tmp) )
4959  {
4960  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4961  WerrorS("Elements of first input ideal must not be equal to -1, 0, 1!");
4962  return TRUE;
4963  }
4964  } else tmp= NULL;
4965  if ( !nIsZero(tmp) )
4966  {
4967  if ( !pIsConstant((p->m)[i]))
4968  {
4969  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4970  WerrorS("Elements of first input ideal must be numbers!");
4971  return TRUE;
4972  }
4973  pevpoint[i]= nCopy( tmp );
4974  }
4975  }
4976 
4977  number *wresults= (number *)omAlloc( m * sizeof( number ) );
4978  for ( i= 0; i < m; i++ )
4979  {
4980  wresults[i]= nInit(0);
4981  if ( (w->m)[i] && !nIsZero(pGetCoeff((w->m)[i])) )
4982  {
4983  if ( !pIsConstant((w->m)[i]))
4984  {
4985  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4986  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4987  WerrorS("Elements of second input ideal must be numbers!");
4988  return TRUE;
4989  }
4990  wresults[i]= nCopy(pGetCoeff((w->m)[i]));
4991  }
4992  }
4993 
4994  vandermonde vm( m, n, tdg, pevpoint, FALSE );
4995  number *ncpoly= vm.interpolateDense( wresults );
4996  // do not free ncpoly[]!!
4997  poly rpoly= vm.numvec2poly( ncpoly );
4998 
4999  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
5000  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
5001 
5002  res->data= (void*)rpoly;
5003  return FALSE;
5004 }
Rational pow(const Rational &a, int e)
Definition: GMPrat.cc:411
int p
Definition: cfModGcd.cc:4080
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