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

4563 {
4564  res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
4565  return FALSE;
4566 }
#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 4568 of file ipshell.cc.

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

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

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

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

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