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

4573 {
4574  res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
4575  return FALSE;
4576 }
#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 4578 of file ipshell.cc.

4579 {
4580  if ( !(rField_is_long_R(currRing)) )
4581  {
4582  WerrorS("Ground field not implemented!");
4583  return TRUE;
4584  }
4585 
4586  simplex * LP;
4587  matrix m;
4588 
4589  leftv v= args;
4590  if ( v->Typ() != MATRIX_CMD ) // 1: matrix
4591  return TRUE;
4592  else
4593  m= (matrix)(v->CopyD());
4594 
4595  LP = new simplex(MATROWS(m),MATCOLS(m));
4596  LP->mapFromMatrix(m);
4597 
4598  v= v->next;
4599  if ( v->Typ() != INT_CMD ) // 2: m = number of constraints
4600  return TRUE;
4601  else
4602  LP->m= (int)(long)(v->Data());
4603 
4604  v= v->next;
4605  if ( v->Typ() != INT_CMD ) // 3: n = number of variables
4606  return TRUE;
4607  else
4608  LP->n= (int)(long)(v->Data());
4609 
4610  v= v->next;
4611  if ( v->Typ() != INT_CMD ) // 4: m1 = number of <= constraints
4612  return TRUE;
4613  else
4614  LP->m1= (int)(long)(v->Data());
4615 
4616  v= v->next;
4617  if ( v->Typ() != INT_CMD ) // 5: m2 = number of >= constraints
4618  return TRUE;
4619  else
4620  LP->m2= (int)(long)(v->Data());
4621 
4622  v= v->next;
4623  if ( v->Typ() != INT_CMD ) // 6: m3 = number of == constraints
4624  return TRUE;
4625  else
4626  LP->m3= (int)(long)(v->Data());
4627 
4628 #ifdef mprDEBUG_PROT
4629  Print("m (constraints) %d\n",LP->m);
4630  Print("n (columns) %d\n",LP->n);
4631  Print("m1 (<=) %d\n",LP->m1);
4632  Print("m2 (>=) %d\n",LP->m2);
4633  Print("m3 (==) %d\n",LP->m3);
4634 #endif
4635 
4636  LP->compute();
4637 
4638  lists lres= (lists)omAlloc( sizeof(slists) );
4639  lres->Init( 6 );
4640 
4641  lres->m[0].rtyp= MATRIX_CMD; // output matrix
4642  lres->m[0].data=(void*)LP->mapToMatrix(m);
4643 
4644  lres->m[1].rtyp= INT_CMD; // found a solution?
4645  lres->m[1].data=(void*)(long)LP->icase;
4646 
4647  lres->m[2].rtyp= INTVEC_CMD;
4648  lres->m[2].data=(void*)LP->posvToIV();
4649 
4650  lres->m[3].rtyp= INTVEC_CMD;
4651  lres->m[3].data=(void*)LP->zrovToIV();
4652 
4653  lres->m[4].rtyp= INT_CMD;
4654  lres->m[4].data=(void*)(long)LP->m;
4655 
4656  lres->m[5].rtyp= INT_CMD;
4657  lres->m[5].data=(void*)(long)LP->n;
4658 
4659  res->data= (void*)lres;
4660 
4661  return FALSE;
4662 }
#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:547
@ 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 4687 of file ipshell.cc.

4688 {
4689  poly gls;
4690  gls= (poly)(arg1->Data());
4691  int howclean= (int)(long)arg3->Data();
4692 
4693  if ( gls == NULL || pIsConstant( gls ) )
4694  {
4695  WerrorS("Input polynomial is constant!");
4696  return TRUE;
4697  }
4698 
4699  if (rField_is_Zp(currRing))
4700  {
4701  int* r=Zp_roots(gls, currRing);
4702  lists rlist;
4703  rlist= (lists)omAlloc( sizeof(slists) );
4704  rlist->Init( r[0] );
4705  for(int i=r[0];i>0;i--)
4706  {
4707  rlist->m[i-1].data=n_Init(r[i],currRing);
4708  rlist->m[i-1].rtyp=NUMBER_CMD;
4709  }
4710  omFree(r);
4711  res->data=rlist;
4712  res->rtyp= LIST_CMD;
4713  return FALSE;
4714  }
4715  if ( !(rField_is_R(currRing) ||
4716  rField_is_Q(currRing) ||
4719  {
4720  WerrorS("Ground field not implemented!");
4721  return TRUE;
4722  }
4723 
4726  {
4727  unsigned long int ii = (unsigned long int)arg2->Data();
4728  setGMPFloatDigits( ii, ii );
4729  }
4730 
4731  int ldummy;
4732  int deg= currRing->pLDeg( gls, &ldummy, currRing );
4733  int i,vpos=0;
4734  poly piter;
4735  lists elist;
4736 
4737  elist= (lists)omAlloc( sizeof(slists) );
4738  elist->Init( 0 );
4739 
4740  if ( rVar(currRing) > 1 )
4741  {
4742  piter= gls;
4743  for ( i= 1; i <= rVar(currRing); i++ )
4744  if ( pGetExp( piter, i ) )
4745  {
4746  vpos= i;
4747  break;
4748  }
4749  while ( piter )
4750  {
4751  for ( i= 1; i <= rVar(currRing); i++ )
4752  if ( (vpos != i) && (pGetExp( piter, i ) != 0) )
4753  {
4754  WerrorS("The input polynomial must be univariate!");
4755  return TRUE;
4756  }
4757  pIter( piter );
4758  }
4759  }
4760 
4761  rootContainer * roots= new rootContainer();
4762  number * pcoeffs= (number *)omAlloc( (deg+1) * sizeof( number ) );
4763  piter= gls;
4764  for ( i= deg; i >= 0; i-- )
4765  {
4766  if ( piter && pTotaldegree(piter) == i )
4767  {
4768  pcoeffs[i]= nCopy( pGetCoeff( piter ) );
4769  //nPrint( pcoeffs[i] );PrintS(" ");
4770  pIter( piter );
4771  }
4772  else
4773  {
4774  pcoeffs[i]= nInit(0);
4775  }
4776  }
4777 
4778 #ifdef mprDEBUG_PROT
4779  for (i=deg; i >= 0; i--)
4780  {
4781  nPrint( pcoeffs[i] );PrintS(" ");
4782  }
4783  PrintLn();
4784 #endif
4785 
4786  roots->fillContainer( pcoeffs, NULL, 1, deg, rootContainer::onepoly, 1 );
4787  roots->solver( howclean );
4788 
4789  int elem= roots->getAnzRoots();
4790  char *dummy;
4791  int j;
4792 
4793  lists rlist;
4794  rlist= (lists)omAlloc( sizeof(slists) );
4795  rlist->Init( elem );
4796 
4798  {
4799  for ( j= 0; j < elem; j++ )
4800  {
4801  rlist->m[j].rtyp=NUMBER_CMD;
4802  rlist->m[j].data=(void *)nCopy((number)(roots->getRoot(j)));
4803  //rlist->m[j].data=(void *)(number)(roots->getRoot(j));
4804  }
4805  }
4806  else
4807  {
4808  for ( j= 0; j < elem; j++ )
4809  {
4810  dummy = complexToStr( (*roots)[j], gmp_output_digits, currRing->cf );
4811  rlist->m[j].rtyp=STRING_CMD;
4812  rlist->m[j].data=(void *)dummy;
4813  }
4814  }
4815 
4816  elist->Clean();
4817  //omFreeSize( (ADDRESS) elist, sizeof(slists) );
4818 
4819  // this is (via fillContainer) the same data as in root
4820  //for ( i= deg; i >= 0; i-- ) nDelete( &pcoeffs[i] );
4821  //omFreeSize( (ADDRESS) pcoeffs, (deg+1) * sizeof( number ) );
4822 
4823  delete roots;
4824 
4825  res->data= (void*)rlist;
4826 
4827  return FALSE;
4828 }
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:523
static BOOLEAN rField_is_Zp(const ring r)
Definition: ring.h:505
static BOOLEAN rField_is_long_C(const ring r)
Definition: ring.h:550
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:511
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:597
@ 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 4664 of file ipshell.cc.

4665 {
4666  ideal gls = (ideal)(arg1->Data());
4667  int imtype= (int)(long)arg2->Data();
4668 
4669  uResultant::resMatType mtype= determineMType( imtype );
4670 
4671  // check input ideal ( = polynomial system )
4672  if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk )
4673  {
4674  return TRUE;
4675  }
4676 
4677  uResultant *resMat= new uResultant( gls, mtype, false );
4678  if (resMat!=NULL)
4679  {
4680  res->rtyp = MODUL_CMD;
4681  res->data= (void*)resMat->accessResMat()->getMatrix();
4682  if (!errorreported) delete resMat;
4683  }
4684  return errorreported;
4685 }
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 4931 of file ipshell.cc.

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

4831 {
4832  int i;
4833  ideal p,w;
4834  p= (ideal)arg1->Data();
4835  w= (ideal)arg2->Data();
4836 
4837  // w[0] = f(p^0)
4838  // w[1] = f(p^1)
4839  // ...
4840  // p can be a vector of numbers (multivariate polynom)
4841  // or one number (univariate polynom)
4842  // tdg = deg(f)
4843 
4844  int n= IDELEMS( p );
4845  int m= IDELEMS( w );
4846  int tdg= (int)(long)arg3->Data();
4847 
4848  res->data= (void*)NULL;
4849 
4850  // check the input
4851  if ( tdg < 1 )
4852  {
4853  WerrorS("Last input parameter must be > 0!");
4854  return TRUE;
4855  }
4856  if ( n != rVar(currRing) )
4857  {
4858  Werror("Size of first input ideal must be equal to %d!",rVar(currRing));
4859  return TRUE;
4860  }
4861  if ( m != (int)pow((double)tdg+1,(double)n) )
4862  {
4863  Werror("Size of second input ideal must be equal to %d!",
4864  (int)pow((double)tdg+1,(double)n));
4865  return TRUE;
4866  }
4867  if ( !(rField_is_Q(currRing) /* ||
4868  rField_is_R() || rField_is_long_R() ||
4869  rField_is_long_C()*/ ) )
4870  {
4871  WerrorS("Ground field not implemented!");
4872  return TRUE;
4873  }
4874 
4875  number tmp;
4876  number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
4877  for ( i= 0; i < n; i++ )
4878  {
4879  pevpoint[i]=nInit(0);
4880  if ( (p->m)[i] )
4881  {
4882  tmp = pGetCoeff( (p->m)[i] );
4883  if ( nIsZero(tmp) || nIsOne(tmp) || nIsMOne(tmp) )
4884  {
4885  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4886  WerrorS("Elements of first input ideal must not be equal to -1, 0, 1!");
4887  return TRUE;
4888  }
4889  } else tmp= NULL;
4890  if ( !nIsZero(tmp) )
4891  {
4892  if ( !pIsConstant((p->m)[i]))
4893  {
4894  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4895  WerrorS("Elements of first input ideal must be numbers!");
4896  return TRUE;
4897  }
4898  pevpoint[i]= nCopy( tmp );
4899  }
4900  }
4901 
4902  number *wresults= (number *)omAlloc( m * sizeof( number ) );
4903  for ( i= 0; i < m; i++ )
4904  {
4905  wresults[i]= nInit(0);
4906  if ( (w->m)[i] && !nIsZero(pGetCoeff((w->m)[i])) )
4907  {
4908  if ( !pIsConstant((w->m)[i]))
4909  {
4910  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4911  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4912  WerrorS("Elements of second input ideal must be numbers!");
4913  return TRUE;
4914  }
4915  wresults[i]= nCopy(pGetCoeff((w->m)[i]));
4916  }
4917  }
4918 
4919  vandermonde vm( m, n, tdg, pevpoint, FALSE );
4920  number *ncpoly= vm.interpolateDense( wresults );
4921  // do not free ncpoly[]!!
4922  poly rpoly= vm.numvec2poly( ncpoly );
4923 
4924  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4925  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4926 
4927  res->data= (void*)rpoly;
4928  return FALSE;
4929 }
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