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

4564 {
4565  res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
4566  return FALSE;
4567 }
#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:3191

◆ loSimplex()

BOOLEAN loSimplex ( leftv  res,
leftv  args 
)

Implementation of the Simplex Algorithm.

For args, see class simplex.

Definition at line 4569 of file ipshell.cc.

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

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

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

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

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