My Project  debian-1:4.1.2-p1+ds-2
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 4537 of file ipshell.cc.

4538 {
4539  res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
4540  return FALSE;
4541 }
#define FALSE
Definition: auxiliary.h:96
void * Data()
Definition: subexpr.cc:1176
CanonicalForm res
Definition: facAbsFact.cc:64
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 4543 of file ipshell.cc.

4544 {
4545  if ( !(rField_is_long_R(currRing)) )
4546  {
4547  WerrorS("Ground field not implemented!");
4548  return TRUE;
4549  }
4550 
4551  simplex * LP;
4552  matrix m;
4553 
4554  leftv v= args;
4555  if ( v->Typ() != MATRIX_CMD ) // 1: matrix
4556  return TRUE;
4557  else
4558  m= (matrix)(v->CopyD());
4559 
4560  LP = new simplex(MATROWS(m),MATCOLS(m));
4561  LP->mapFromMatrix(m);
4562 
4563  v= v->next;
4564  if ( v->Typ() != INT_CMD ) // 2: m = number of constraints
4565  return TRUE;
4566  else
4567  LP->m= (int)(long)(v->Data());
4568 
4569  v= v->next;
4570  if ( v->Typ() != INT_CMD ) // 3: n = number of variables
4571  return TRUE;
4572  else
4573  LP->n= (int)(long)(v->Data());
4574 
4575  v= v->next;
4576  if ( v->Typ() != INT_CMD ) // 4: m1 = number of <= constraints
4577  return TRUE;
4578  else
4579  LP->m1= (int)(long)(v->Data());
4580 
4581  v= v->next;
4582  if ( v->Typ() != INT_CMD ) // 5: m2 = number of >= constraints
4583  return TRUE;
4584  else
4585  LP->m2= (int)(long)(v->Data());
4586 
4587  v= v->next;
4588  if ( v->Typ() != INT_CMD ) // 6: m3 = number of == constraints
4589  return TRUE;
4590  else
4591  LP->m3= (int)(long)(v->Data());
4592 
4593 #ifdef mprDEBUG_PROT
4594  Print("m (constraints) %d\n",LP->m);
4595  Print("n (columns) %d\n",LP->n);
4596  Print("m1 (<=) %d\n",LP->m1);
4597  Print("m2 (>=) %d\n",LP->m2);
4598  Print("m3 (==) %d\n",LP->m3);
4599 #endif
4600 
4601  LP->compute();
4602 
4603  lists lres= (lists)omAlloc( sizeof(slists) );
4604  lres->Init( 6 );
4605 
4606  lres->m[0].rtyp= MATRIX_CMD; // output matrix
4607  lres->m[0].data=(void*)LP->mapToMatrix(m);
4608 
4609  lres->m[1].rtyp= INT_CMD; // found a solution?
4610  lres->m[1].data=(void*)(long)LP->icase;
4611 
4612  lres->m[2].rtyp= INTVEC_CMD;
4613  lres->m[2].data=(void*)LP->posvToIV();
4614 
4615  lres->m[3].rtyp= INTVEC_CMD;
4616  lres->m[3].data=(void*)LP->zrovToIV();
4617 
4618  lres->m[4].rtyp= INT_CMD;
4619  lres->m[4].data=(void*)(long)LP->m;
4620 
4621  lres->m[5].rtyp= INT_CMD;
4622  lres->m[5].data=(void*)(long)LP->n;
4623 
4624  res->data= (void*)lres;
4625 
4626  return FALSE;
4627 }
#define TRUE
Definition: auxiliary.h:100
int m
Definition: cfEzgcd.cc:121
Variable next() const
Definition: factory.h:137
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:37
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:537
@ 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 4652 of file ipshell.cc.

4653 {
4654 
4655  poly gls;
4656  gls= (poly)(arg1->Data());
4657  int howclean= (int)(long)arg3->Data();
4658 
4659  if ( !(rField_is_R(currRing) ||
4660  rField_is_Q(currRing) ||
4663  {
4664  WerrorS("Ground field not implemented!");
4665  return TRUE;
4666  }
4667 
4670  {
4671  unsigned long int ii = (unsigned long int)arg2->Data();
4672  setGMPFloatDigits( ii, ii );
4673  }
4674 
4675  if ( gls == NULL || pIsConstant( gls ) )
4676  {
4677  WerrorS("Input polynomial is constant!");
4678  return TRUE;
4679  }
4680 
4681  int ldummy;
4682  int deg= currRing->pLDeg( gls, &ldummy, currRing );
4683  int i,vpos=0;
4684  poly piter;
4685  lists elist;
4686  lists rlist;
4687 
4688  elist= (lists)omAlloc( sizeof(slists) );
4689  elist->Init( 0 );
4690 
4691  if ( rVar(currRing) > 1 )
4692  {
4693  piter= gls;
4694  for ( i= 1; i <= rVar(currRing); i++ )
4695  if ( pGetExp( piter, i ) )
4696  {
4697  vpos= i;
4698  break;
4699  }
4700  while ( piter )
4701  {
4702  for ( i= 1; i <= rVar(currRing); i++ )
4703  if ( (vpos != i) && (pGetExp( piter, i ) != 0) )
4704  {
4705  WerrorS("The input polynomial must be univariate!");
4706  return TRUE;
4707  }
4708  pIter( piter );
4709  }
4710  }
4711 
4712  rootContainer * roots= new rootContainer();
4713  number * pcoeffs= (number *)omAlloc( (deg+1) * sizeof( number ) );
4714  piter= gls;
4715  for ( i= deg; i >= 0; i-- )
4716  {
4717  if ( piter && pTotaldegree(piter) == i )
4718  {
4719  pcoeffs[i]= nCopy( pGetCoeff( piter ) );
4720  //nPrint( pcoeffs[i] );PrintS(" ");
4721  pIter( piter );
4722  }
4723  else
4724  {
4725  pcoeffs[i]= nInit(0);
4726  }
4727  }
4728 
4729 #ifdef mprDEBUG_PROT
4730  for (i=deg; i >= 0; i--)
4731  {
4732  nPrint( pcoeffs[i] );PrintS(" ");
4733  }
4734  PrintLn();
4735 #endif
4736 
4737  roots->fillContainer( pcoeffs, NULL, 1, deg, rootContainer::onepoly, 1 );
4738  roots->solver( howclean );
4739 
4740  int elem= roots->getAnzRoots();
4741  char *dummy;
4742  int j;
4743 
4744  rlist= (lists)omAlloc( sizeof(slists) );
4745  rlist->Init( elem );
4746 
4748  {
4749  for ( j= 0; j < elem; j++ )
4750  {
4751  rlist->m[j].rtyp=NUMBER_CMD;
4752  rlist->m[j].data=(void *)nCopy((number)(roots->getRoot(j)));
4753  //rlist->m[j].data=(void *)(number)(roots->getRoot(j));
4754  }
4755  }
4756  else
4757  {
4758  for ( j= 0; j < elem; j++ )
4759  {
4760  dummy = complexToStr( (*roots)[j], gmp_output_digits, currRing->cf );
4761  rlist->m[j].rtyp=STRING_CMD;
4762  rlist->m[j].data=(void *)dummy;
4763  }
4764  }
4765 
4766  elist->Clean();
4767  //omFreeSize( (ADDRESS) elist, sizeof(slists) );
4768 
4769  // this is (via fillContainer) the same data as in root
4770  //for ( i= deg; i >= 0; i-- ) nDelete( &pcoeffs[i] );
4771  //omFreeSize( (ADDRESS) pcoeffs, (deg+1) * sizeof( number ) );
4772 
4773  delete roots;
4774 
4775  res->rtyp= LIST_CMD;
4776  res->data= (void*)rlist;
4777 
4778  return FALSE;
4779 }
int i
Definition: cfEzgcd.cc:125
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
int j
Definition: facHensel.cc:105
@ 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 NULL
Definition: omList.c:12
static long pTotaldegree(poly p)
Definition: polys.h:278
#define pIsConstant(p)
like above, except that Comp must be 0
Definition: polys.h:234
#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:513
static BOOLEAN rField_is_long_C(const ring r)
Definition: ring.h:540
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:501
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:587
@ LIST_CMD
Definition: tok.h:118
@ STRING_CMD
Definition: tok.h:184

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

4630 {
4631  ideal gls = (ideal)(arg1->Data());
4632  int imtype= (int)(long)arg2->Data();
4633 
4634  uResultant::resMatType mtype= determineMType( imtype );
4635 
4636  // check input ideal ( = polynomial system )
4637  if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk )
4638  {
4639  return TRUE;
4640  }
4641 
4642  uResultant *resMat= new uResultant( gls, mtype, false );
4643  if (resMat!=NULL)
4644  {
4645  res->rtyp = MODUL_CMD;
4646  res->data= (void*)resMat->accessResMat()->getMatrix();
4647  if (!errorreported) delete resMat;
4648  }
4649  return errorreported;
4650 }
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 4882 of file ipshell.cc.

4883 {
4884  leftv v= args;
4885 
4886  ideal gls;
4887  int imtype;
4888  int howclean;
4889 
4890  // get ideal
4891  if ( v->Typ() != IDEAL_CMD )
4892  return TRUE;
4893  else gls= (ideal)(v->Data());
4894  v= v->next;
4895 
4896  // get resultant matrix type to use (0,1)
4897  if ( v->Typ() != INT_CMD )
4898  return TRUE;
4899  else imtype= (int)(long)v->Data();
4900  v= v->next;
4901 
4902  if (imtype==0)
4903  {
4904  ideal test_id=idInit(1,1);
4905  int j;
4906  for(j=IDELEMS(gls)-1;j>=0;j--)
4907  {
4908  if (gls->m[j]!=NULL)
4909  {
4910  test_id->m[0]=gls->m[j];
4911  intvec *dummy_w=id_QHomWeight(test_id, currRing);
4912  if (dummy_w!=NULL)
4913  {
4914  WerrorS("Newton polytope not of expected dimension");
4915  delete dummy_w;
4916  return TRUE;
4917  }
4918  }
4919  }
4920  }
4921 
4922  // get and set precision in digits ( > 0 )
4923  if ( v->Typ() != INT_CMD )
4924  return TRUE;
4925  else if ( !(rField_is_R(currRing) || rField_is_long_R(currRing) || \
4927  {
4928  unsigned long int ii=(unsigned long int)v->Data();
4929  setGMPFloatDigits( ii, ii );
4930  }
4931  v= v->next;
4932 
4933  // get interpolation steps (0,1,2)
4934  if ( v->Typ() != INT_CMD )
4935  return TRUE;
4936  else howclean= (int)(long)v->Data();
4937 
4938  uResultant::resMatType mtype= determineMType( imtype );
4939  int i,count;
4940  lists listofroots= NULL;
4941  number smv= NULL;
4942  BOOLEAN interpolate_det= (mtype==uResultant::denseResMat)?TRUE:FALSE;
4943 
4944  //emptylist= (lists)omAlloc( sizeof(slists) );
4945  //emptylist->Init( 0 );
4946 
4947  //res->rtyp = LIST_CMD;
4948  //res->data= (void *)emptylist;
4949 
4950  // check input ideal ( = polynomial system )
4951  if ( mprIdealCheck( gls, args->Name(), mtype ) != mprOk )
4952  {
4953  return TRUE;
4954  }
4955 
4956  uResultant * ures;
4957  rootContainer ** iproots;
4958  rootContainer ** muiproots;
4959  rootArranger * arranger;
4960 
4961  // main task 1: setup of resultant matrix
4962  ures= new uResultant( gls, mtype );
4963  if ( ures->accessResMat()->initState() != resMatrixBase::ready )
4964  {
4965  WerrorS("Error occurred during matrix setup!");
4966  return TRUE;
4967  }
4968 
4969  // if dense resultant, check if minor nonsingular
4970  if ( mtype == uResultant::denseResMat )
4971  {
4972  smv= ures->accessResMat()->getSubDet();
4973 #ifdef mprDEBUG_PROT
4974  PrintS("// Determinant of submatrix: ");nPrint(smv);PrintLn();
4975 #endif
4976  if ( nIsZero(smv) )
4977  {
4978  WerrorS("Unsuitable input ideal: Minor of resultant matrix is singular!");
4979  return TRUE;
4980  }
4981  }
4982 
4983  // main task 2: Interpolate specialized resultant polynomials
4984  if ( interpolate_det )
4985  iproots= ures->interpolateDenseSP( false, smv );
4986  else
4987  iproots= ures->specializeInU( false, smv );
4988 
4989  // main task 3: Interpolate specialized resultant polynomials
4990  if ( interpolate_det )
4991  muiproots= ures->interpolateDenseSP( true, smv );
4992  else
4993  muiproots= ures->specializeInU( true, smv );
4994 
4995 #ifdef mprDEBUG_PROT
4996  int c= iproots[0]->getAnzElems();
4997  for (i=0; i < c; i++) pWrite(iproots[i]->getPoly());
4998  c= muiproots[0]->getAnzElems();
4999  for (i=0; i < c; i++) pWrite(muiproots[i]->getPoly());
5000 #endif
5001 
5002  // main task 4: Compute roots of specialized polys and match them up
5003  arranger= new rootArranger( iproots, muiproots, howclean );
5004  arranger->solve_all();
5005 
5006  // get list of roots
5007  if ( arranger->success() )
5008  {
5009  arranger->arrange();
5010  listofroots= listOfRoots(arranger, gmp_output_digits );
5011  }
5012  else
5013  {
5014  WerrorS("Solver was unable to find any roots!");
5015  return TRUE;
5016  }
5017 
5018  // free everything
5019  count= iproots[0]->getAnzElems();
5020  for (i=0; i < count; i++) delete iproots[i];
5021  omFreeSize( (ADDRESS) iproots, count * sizeof(rootContainer*) );
5022  count= muiproots[0]->getAnzElems();
5023  for (i=0; i < count; i++) delete muiproots[i];
5024  omFreeSize( (ADDRESS) muiproots, count * sizeof(rootContainer*) );
5025 
5026  delete ures;
5027  delete arranger;
5028  nDelete( &smv );
5029 
5030  res->data= (void *)listofroots;
5031 
5032  //emptylist->Clean();
5033  // omFreeSize( (ADDRESS) emptylist, sizeof(slists) );
5034 
5035  return FALSE;
5036 }
int BOOLEAN
Definition: auxiliary.h:87
void * ADDRESS
Definition: auxiliary.h:135
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:5039
#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:304
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 4781 of file ipshell.cc.

4782 {
4783  int i;
4784  ideal p,w;
4785  p= (ideal)arg1->Data();
4786  w= (ideal)arg2->Data();
4787 
4788  // w[0] = f(p^0)
4789  // w[1] = f(p^1)
4790  // ...
4791  // p can be a vector of numbers (multivariate polynom)
4792  // or one number (univariate polynom)
4793  // tdg = deg(f)
4794 
4795  int n= IDELEMS( p );
4796  int m= IDELEMS( w );
4797  int tdg= (int)(long)arg3->Data();
4798 
4799  res->data= (void*)NULL;
4800 
4801  // check the input
4802  if ( tdg < 1 )
4803  {
4804  WerrorS("Last input parameter must be > 0!");
4805  return TRUE;
4806  }
4807  if ( n != rVar(currRing) )
4808  {
4809  Werror("Size of first input ideal must be equal to %d!",rVar(currRing));
4810  return TRUE;
4811  }
4812  if ( m != (int)pow((double)tdg+1,(double)n) )
4813  {
4814  Werror("Size of second input ideal must be equal to %d!",
4815  (int)pow((double)tdg+1,(double)n));
4816  return TRUE;
4817  }
4818  if ( !(rField_is_Q(currRing) /* ||
4819  rField_is_R() || rField_is_long_R() ||
4820  rField_is_long_C()*/ ) )
4821  {
4822  WerrorS("Ground field not implemented!");
4823  return TRUE;
4824  }
4825 
4826  number tmp;
4827  number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
4828  for ( i= 0; i < n; i++ )
4829  {
4830  pevpoint[i]=nInit(0);
4831  if ( (p->m)[i] )
4832  {
4833  tmp = pGetCoeff( (p->m)[i] );
4834  if ( nIsZero(tmp) || nIsOne(tmp) || nIsMOne(tmp) )
4835  {
4836  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4837  WerrorS("Elements of first input ideal must not be equal to -1, 0, 1!");
4838  return TRUE;
4839  }
4840  } else tmp= NULL;
4841  if ( !nIsZero(tmp) )
4842  {
4843  if ( !pIsConstant((p->m)[i]))
4844  {
4845  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4846  WerrorS("Elements of first input ideal must be numbers!");
4847  return TRUE;
4848  }
4849  pevpoint[i]= nCopy( tmp );
4850  }
4851  }
4852 
4853  number *wresults= (number *)omAlloc( m * sizeof( number ) );
4854  for ( i= 0; i < m; i++ )
4855  {
4856  wresults[i]= nInit(0);
4857  if ( (w->m)[i] && !nIsZero(pGetCoeff((w->m)[i])) )
4858  {
4859  if ( !pIsConstant((w->m)[i]))
4860  {
4861  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4862  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4863  WerrorS("Elements of second input ideal must be numbers!");
4864  return TRUE;
4865  }
4866  wresults[i]= nCopy(pGetCoeff((w->m)[i]));
4867  }
4868  }
4869 
4870  vandermonde vm( m, n, tdg, pevpoint, FALSE );
4871  number *ncpoly= vm.interpolateDense( wresults );
4872  // do not free ncpoly[]!!
4873  poly rpoly= vm.numvec2poly( ncpoly );
4874 
4875  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4876  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4877 
4878  res->data= (void*)rpoly;
4879  return FALSE;
4880 }
Rational pow(const Rational &a, int e)
Definition: GMPrat.cc:411
int p
Definition: cfModGcd.cc:4019
vandermonde system solver for interpolating polynomials from their values
Definition: mpr_numeric.h:29
const CanonicalForm & w
Definition: facAbsFact.cc:55
#define nIsMOne(n)
Definition: numbers.h:26
#define nIsOne(n)
Definition: numbers.h:25
void Werror(const char *fmt,...)
Definition: reporter.cc:189