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

4567 {
4568  res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
4569  return FALSE;
4570 }
#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 4572 of file ipshell.cc.

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

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

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

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

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