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

#define DEFAULT_DIGITS   30

Definition at line 13 of file mpr_inout.h.

#define MPR_DENSE   1

Definition at line 15 of file mpr_inout.h.

#define MPR_SPARSE   2

Definition at line 16 of file mpr_inout.h.

Function Documentation

BOOLEAN loNewtonP ( leftv  res,
leftv  arg1 
)

compute Newton Polytopes of input polynomials

Definition at line 4478 of file ipshell.cc.

4479 {
4480  res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
4481  return FALSE;
4482 }
#define FALSE
Definition: auxiliary.h:140
ideal loNewtonPolytope(const ideal id)
Definition: mpr_base.cc:3192
void * data
Definition: subexpr.h:89
void * Data()
Definition: subexpr.cc:1118
BOOLEAN loSimplex ( leftv  res,
leftv  args 
)

Implementation of the Simplex Algorithm.

For args, see class simplex.

Definition at line 4484 of file ipshell.cc.

4485 {
4486  if ( !(rField_is_long_R(currRing)) )
4487  {
4488  WerrorS("Ground field not implemented!");
4489  return TRUE;
4490  }
4491 
4492  simplex * LP;
4493  matrix m;
4494 
4495  leftv v= args;
4496  if ( v->Typ() != MATRIX_CMD ) // 1: matrix
4497  return TRUE;
4498  else
4499  m= (matrix)(v->CopyD());
4500 
4501  LP = new simplex(MATROWS(m),MATCOLS(m));
4502  LP->mapFromMatrix(m);
4503 
4504  v= v->next;
4505  if ( v->Typ() != INT_CMD ) // 2: m = number of constraints
4506  return TRUE;
4507  else
4508  LP->m= (int)(long)(v->Data());
4509 
4510  v= v->next;
4511  if ( v->Typ() != INT_CMD ) // 3: n = number of variables
4512  return TRUE;
4513  else
4514  LP->n= (int)(long)(v->Data());
4515 
4516  v= v->next;
4517  if ( v->Typ() != INT_CMD ) // 4: m1 = number of <= constraints
4518  return TRUE;
4519  else
4520  LP->m1= (int)(long)(v->Data());
4521 
4522  v= v->next;
4523  if ( v->Typ() != INT_CMD ) // 5: m2 = number of >= constraints
4524  return TRUE;
4525  else
4526  LP->m2= (int)(long)(v->Data());
4527 
4528  v= v->next;
4529  if ( v->Typ() != INT_CMD ) // 6: m3 = number of == constraints
4530  return TRUE;
4531  else
4532  LP->m3= (int)(long)(v->Data());
4533 
4534 #ifdef mprDEBUG_PROT
4535  Print("m (constraints) %d\n",LP->m);
4536  Print("n (columns) %d\n",LP->n);
4537  Print("m1 (<=) %d\n",LP->m1);
4538  Print("m2 (>=) %d\n",LP->m2);
4539  Print("m3 (==) %d\n",LP->m3);
4540 #endif
4541 
4542  LP->compute();
4543 
4544  lists lres= (lists)omAlloc( sizeof(slists) );
4545  lres->Init( 6 );
4546 
4547  lres->m[0].rtyp= MATRIX_CMD; // output matrix
4548  lres->m[0].data=(void*)LP->mapToMatrix(m);
4549 
4550  lres->m[1].rtyp= INT_CMD; // found a solution?
4551  lres->m[1].data=(void*)(long)LP->icase;
4552 
4553  lres->m[2].rtyp= INTVEC_CMD;
4554  lres->m[2].data=(void*)LP->posvToIV();
4555 
4556  lres->m[3].rtyp= INTVEC_CMD;
4557  lres->m[3].data=(void*)LP->zrovToIV();
4558 
4559  lres->m[4].rtyp= INT_CMD;
4560  lres->m[4].data=(void*)(long)LP->m;
4561 
4562  lres->m[5].rtyp= INT_CMD;
4563  lres->m[5].data=(void*)(long)LP->n;
4564 
4565  res->data= (void*)lres;
4566 
4567  return FALSE;
4568 }
sleftv * m
Definition: lists.h:45
Class used for (list of) interpreter objects.
Definition: subexpr.h:83
matrix mapToMatrix(matrix m)
void compute()
#define Print
Definition: emacs.cc:83
Definition: tok.h:98
Definition: lists.h:22
#define FALSE
Definition: auxiliary.h:140
Linear Programming / Linear Optimization using Simplex - Algorithm.
Definition: mpr_numeric.h:194
#define TRUE
Definition: auxiliary.h:144
intvec * zrovToIV()
void WerrorS(const char *s)
Definition: feFopen.cc:24
int Typ()
Definition: subexpr.cc:976
#define omAlloc(size)
Definition: omAllocDecl.h:210
void * data
Definition: subexpr.h:89
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:12
intvec * posvToIV()
BOOLEAN mapFromMatrix(matrix m)
ip_smatrix * matrix
int m
Definition: cfEzgcd.cc:119
leftv next
Definition: subexpr.h:87
INLINE_THIS void Init(int l=0)
Definition: lists.h:66
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
#define MATCOLS(i)
Definition: matpol.h:28
slists * lists
Definition: mpr_numeric.h:146
static BOOLEAN rField_is_long_R(const ring r)
Definition: ring.h:488
int rtyp
Definition: subexpr.h:92
void * Data()
Definition: subexpr.cc:1118
#define MATROWS(i)
Definition: matpol.h:27
int icase
Definition: mpr_numeric.h:201
void * CopyD(int t)
Definition: subexpr.cc:676
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 4593 of file ipshell.cc.

4594 {
4595 
4596  poly gls;
4597  gls= (poly)(arg1->Data());
4598  int howclean= (int)(long)arg3->Data();
4599 
4600  if ( !(rField_is_R(currRing) ||
4601  rField_is_Q(currRing) ||
4604  {
4605  WerrorS("Ground field not implemented!");
4606  return TRUE;
4607  }
4608 
4611  {
4612  unsigned long int ii = (unsigned long int)arg2->Data();
4613  setGMPFloatDigits( ii, ii );
4614  }
4615 
4616  if ( gls == NULL || pIsConstant( gls ) )
4617  {
4618  WerrorS("Input polynomial is constant!");
4619  return TRUE;
4620  }
4621 
4622  int ldummy;
4623  int deg= currRing->pLDeg( gls, &ldummy, currRing );
4624  // int deg= pDeg( gls );
4625  // int len= pLength( gls );
4626  int i,vpos=0;
4627  poly piter;
4628  lists elist;
4629  lists rlist;
4630 
4631  elist= (lists)omAlloc( sizeof(slists) );
4632  elist->Init( 0 );
4633 
4634  if ( rVar(currRing) > 1 )
4635  {
4636  piter= gls;
4637  for ( i= 1; i <= rVar(currRing); i++ )
4638  if ( pGetExp( piter, i ) )
4639  {
4640  vpos= i;
4641  break;
4642  }
4643  while ( piter )
4644  {
4645  for ( i= 1; i <= rVar(currRing); i++ )
4646  if ( (vpos != i) && (pGetExp( piter, i ) != 0) )
4647  {
4648  WerrorS("The input polynomial must be univariate!");
4649  return TRUE;
4650  }
4651  pIter( piter );
4652  }
4653  }
4654 
4655  rootContainer * roots= new rootContainer();
4656  number * pcoeffs= (number *)omAlloc( (deg+1) * sizeof( number ) );
4657  piter= gls;
4658  for ( i= deg; i >= 0; i-- )
4659  {
4660  //if ( piter ) Print("deg %d, pDeg(piter) %d\n",i,pTotaldegree(piter));
4661  if ( piter && pTotaldegree(piter) == i )
4662  {
4663  pcoeffs[i]= nCopy( pGetCoeff( piter ) );
4664  //nPrint( pcoeffs[i] );PrintS(" ");
4665  pIter( piter );
4666  }
4667  else
4668  {
4669  pcoeffs[i]= nInit(0);
4670  }
4671  }
4672 
4673 #ifdef mprDEBUG_PROT
4674  for (i=deg; i >= 0; i--)
4675  {
4676  nPrint( pcoeffs[i] );PrintS(" ");
4677  }
4678  PrintLn();
4679 #endif
4680 
4681  roots->fillContainer( pcoeffs, NULL, 1, deg, rootContainer::onepoly, 1 );
4682  roots->solver( howclean );
4683 
4684  int elem= roots->getAnzRoots();
4685  char *dummy;
4686  int j;
4687 
4688  rlist= (lists)omAlloc( sizeof(slists) );
4689  rlist->Init( elem );
4690 
4692  {
4693  for ( j= 0; j < elem; j++ )
4694  {
4695  rlist->m[j].rtyp=NUMBER_CMD;
4696  rlist->m[j].data=(void *)nCopy((number)(roots->getRoot(j)));
4697  //rlist->m[j].data=(void *)(number)(roots->getRoot(j));
4698  }
4699  }
4700  else
4701  {
4702  for ( j= 0; j < elem; j++ )
4703  {
4704  dummy = complexToStr( (*roots)[j], gmp_output_digits, currRing->cf );
4705  rlist->m[j].rtyp=STRING_CMD;
4706  rlist->m[j].data=(void *)dummy;
4707  }
4708  }
4709 
4710  elist->Clean();
4711  //omFreeSize( (ADDRESS) elist, sizeof(slists) );
4712 
4713  // this is (via fillContainer) the same data as in root
4714  //for ( i= deg; i >= 0; i-- ) nDelete( &pcoeffs[i] );
4715  //omFreeSize( (ADDRESS) pcoeffs, (deg+1) * sizeof( number ) );
4716 
4717  delete roots;
4718 
4719  res->rtyp= LIST_CMD;
4720  res->data= (void*)rlist;
4721 
4722  return FALSE;
4723 }
complex root finder for univariate polynomials based on laguers algorithm
Definition: mpr_numeric.h:65
sleftv * m
Definition: lists.h:45
void PrintLn()
Definition: reporter.cc:327
Definition: lists.h:22
#define FALSE
Definition: auxiliary.h:140
static BOOLEAN rField_is_R(const ring r)
Definition: ring.h:464
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:537
#define TRUE
Definition: auxiliary.h:144
bool solver(const int polishmode=PM_NONE)
Definition: mpr_numeric.cc:450
void WerrorS(const char *s)
Definition: feFopen.cc:24
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy ...
Definition: monomials.h:51
#define omAlloc(size)
Definition: omAllocDecl.h:210
#define nPrint(a)
only for debug, over any initalized currRing
Definition: numbers.h:46
void * data
Definition: subexpr.h:89
#define pIter(p)
Definition: monomials.h:44
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:12
#define pGetExp(p, i)
Exponent.
Definition: polys.h:41
int j
Definition: myNF.cc:70
static long pTotaldegree(poly p)
Definition: polys.h:253
#define pIsConstant(p)
like above, except that Comp might be != 0
Definition: polys.h:209
void fillContainer(number *_coeffs, number *_ievpoint, const int _var, const int _tdg, const rootType _rt, const int _anz)
Definition: mpr_numeric.cc:313
int i
Definition: cfEzgcd.cc:123
void PrintS(const char *s)
Definition: reporter.cc:294
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:458
gmp_complex * getRoot(const int i)
Definition: mpr_numeric.h:88
static BOOLEAN rField_is_long_C(const ring r)
Definition: ring.h:491
INLINE_THIS void Init(int l=0)
Definition: lists.h:66
#define NULL
Definition: omList.c:10
slists * lists
Definition: mpr_numeric.h:146
int getAnzRoots()
Definition: mpr_numeric.h:97
static BOOLEAN rField_is_long_R(const ring r)
Definition: ring.h:488
int rtyp
Definition: subexpr.h:92
#define nCopy(n)
Definition: numbers.h:15
void Clean(ring r=currRing)
Definition: lists.h:25
void * Data()
Definition: subexpr.cc:1118
Definition: tok.h:120
char * complexToStr(gmp_complex &c, const unsigned int oprec, const coeffs src)
Definition: mpr_complex.cc:717
size_t gmp_output_digits
Definition: mpr_complex.cc:44
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:62
polyrec * poly
Definition: hilb.h:10
#define nInit(i)
Definition: numbers.h:24
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 4570 of file ipshell.cc.

4571 {
4572  ideal gls = (ideal)(arg1->Data());
4573  int imtype= (int)(long)arg2->Data();
4574 
4575  uResultant::resMatType mtype= determineMType( imtype );
4576 
4577  // check input ideal ( = polynomial system )
4578  if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk )
4579  {
4580  return TRUE;
4581  }
4582 
4583  uResultant *resMat= new uResultant( gls, mtype, false );
4584  if (resMat!=NULL)
4585  {
4586  res->rtyp = MODUL_CMD;
4587  res->data= (void*)resMat->accessResMat()->getMatrix();
4588  if (!errorreported) delete resMat;
4589  }
4590  return errorreported;
4591 }
Base class for solving 0-dim poly systems using u-resultant.
Definition: mpr_base.h:62
resMatrixBase * accessResMat()
Definition: mpr_base.h:78
#define TRUE
Definition: auxiliary.h:144
uResultant::resMatType determineMType(int imtype)
const char * Name()
Definition: subexpr.h:121
Definition: mpr_base.h:98
void * data
Definition: subexpr.h:89
virtual ideal getMatrix()
Definition: mpr_base.h:31
mprState mprIdealCheck(const ideal theIdeal, const char *name, uResultant::resMatType mtype, BOOLEAN rmatrix=false)
short errorreported
Definition: feFopen.cc:23
#define NULL
Definition: omList.c:10
int rtyp
Definition: subexpr.h:92
void * Data()
Definition: subexpr.cc:1118
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 4826 of file ipshell.cc.

4827 {
4828  leftv v= args;
4829 
4830  ideal gls;
4831  int imtype;
4832  int howclean;
4833 
4834  // get ideal
4835  if ( v->Typ() != IDEAL_CMD )
4836  return TRUE;
4837  else gls= (ideal)(v->Data());
4838  v= v->next;
4839 
4840  // get resultant matrix type to use (0,1)
4841  if ( v->Typ() != INT_CMD )
4842  return TRUE;
4843  else imtype= (int)(long)v->Data();
4844  v= v->next;
4845 
4846  if (imtype==0)
4847  {
4848  ideal test_id=idInit(1,1);
4849  int j;
4850  for(j=IDELEMS(gls)-1;j>=0;j--)
4851  {
4852  if (gls->m[j]!=NULL)
4853  {
4854  test_id->m[0]=gls->m[j];
4855  intvec *dummy_w=id_QHomWeight(test_id, currRing);
4856  if (dummy_w!=NULL)
4857  {
4858  WerrorS("Newton polytope not of expected dimension");
4859  delete dummy_w;
4860  return TRUE;
4861  }
4862  }
4863  }
4864  }
4865 
4866  // get and set precision in digits ( > 0 )
4867  if ( v->Typ() != INT_CMD )
4868  return TRUE;
4869  else if ( !(rField_is_R(currRing) || rField_is_long_R(currRing) || \
4871  {
4872  unsigned long int ii=(unsigned long int)v->Data();
4873  setGMPFloatDigits( ii, ii );
4874  }
4875  v= v->next;
4876 
4877  // get interpolation steps (0,1,2)
4878  if ( v->Typ() != INT_CMD )
4879  return TRUE;
4880  else howclean= (int)(long)v->Data();
4881 
4882  uResultant::resMatType mtype= determineMType( imtype );
4883  int i,count;
4884  lists listofroots= NULL;
4885  number smv= NULL;
4886  BOOLEAN interpolate_det= (mtype==uResultant::denseResMat)?TRUE:FALSE;
4887 
4888  //emptylist= (lists)omAlloc( sizeof(slists) );
4889  //emptylist->Init( 0 );
4890 
4891  //res->rtyp = LIST_CMD;
4892  //res->data= (void *)emptylist;
4893 
4894  // check input ideal ( = polynomial system )
4895  if ( mprIdealCheck( gls, args->Name(), mtype ) != mprOk )
4896  {
4897  return TRUE;
4898  }
4899 
4900  uResultant * ures;
4901  rootContainer ** iproots;
4902  rootContainer ** muiproots;
4903  rootArranger * arranger;
4904 
4905  // main task 1: setup of resultant matrix
4906  ures= new uResultant( gls, mtype );
4907  if ( ures->accessResMat()->initState() != resMatrixBase::ready )
4908  {
4909  WerrorS("Error occurred during matrix setup!");
4910  return TRUE;
4911  }
4912 
4913  // if dense resultant, check if minor nonsingular
4914  if ( mtype == uResultant::denseResMat )
4915  {
4916  smv= ures->accessResMat()->getSubDet();
4917 #ifdef mprDEBUG_PROT
4918  PrintS("// Determinant of submatrix: ");nPrint(smv);PrintLn();
4919 #endif
4920  if ( nIsZero(smv) )
4921  {
4922  WerrorS("Unsuitable input ideal: Minor of resultant matrix is singular!");
4923  return TRUE;
4924  }
4925  }
4926 
4927  // main task 2: Interpolate specialized resultant polynomials
4928  if ( interpolate_det )
4929  iproots= ures->interpolateDenseSP( false, smv );
4930  else
4931  iproots= ures->specializeInU( false, smv );
4932 
4933  // main task 3: Interpolate specialized resultant polynomials
4934  if ( interpolate_det )
4935  muiproots= ures->interpolateDenseSP( true, smv );
4936  else
4937  muiproots= ures->specializeInU( true, smv );
4938 
4939 #ifdef mprDEBUG_PROT
4940  int c= iproots[0]->getAnzElems();
4941  for (i=0; i < c; i++) pWrite(iproots[i]->getPoly());
4942  c= muiproots[0]->getAnzElems();
4943  for (i=0; i < c; i++) pWrite(muiproots[i]->getPoly());
4944 #endif
4945 
4946  // main task 4: Compute roots of specialized polys and match them up
4947  arranger= new rootArranger( iproots, muiproots, howclean );
4948  arranger->solve_all();
4949 
4950  // get list of roots
4951  if ( arranger->success() )
4952  {
4953  arranger->arrange();
4954  listofroots= listOfRoots(arranger, gmp_output_digits );
4955  }
4956  else
4957  {
4958  WerrorS("Solver was unable to find any roots!");
4959  return TRUE;
4960  }
4961 
4962  // free everything
4963  count= iproots[0]->getAnzElems();
4964  for (i=0; i < count; i++) delete iproots[i];
4965  omFreeSize( (ADDRESS) iproots, count * sizeof(rootContainer*) );
4966  count= muiproots[0]->getAnzElems();
4967  for (i=0; i < count; i++) delete muiproots[i];
4968  omFreeSize( (ADDRESS) muiproots, count * sizeof(rootContainer*) );
4969 
4970  delete ures;
4971  delete arranger;
4972  nDelete( &smv );
4973 
4974  res->data= (void *)listofroots;
4975 
4976  //emptylist->Clean();
4977  // omFreeSize( (ADDRESS) emptylist, sizeof(slists) );
4978 
4979  return FALSE;
4980 }
int status int void size_t count
Definition: si_signals.h:59
complex root finder for univariate polynomials based on laguers algorithm
Definition: mpr_numeric.h:65
Class used for (list of) interpreter objects.
Definition: subexpr.h:83
void PrintLn()
Definition: reporter.cc:327
Base class for solving 0-dim poly systems using u-resultant.
Definition: mpr_base.h:62
Definition: tok.h:98
Definition: lists.h:22
virtual IStateType initState() const
Definition: mpr_base.h:41
#define FALSE
Definition: auxiliary.h:140
static BOOLEAN rField_is_R(const ring r)
Definition: ring.h:464
resMatrixBase * accessResMat()
Definition: mpr_base.h:78
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
intvec * id_QHomWeight(ideal id, const ring r)
#define TRUE
Definition: auxiliary.h:144
uResultant::resMatType determineMType(int imtype)
void * ADDRESS
Definition: auxiliary.h:161
void pWrite(poly p)
Definition: polys.h:279
void WerrorS(const char *s)
Definition: feFopen.cc:24
int getAnzElems()
Definition: mpr_numeric.h:95
rootContainer ** specializeInU(BOOLEAN matchUp=false, const number subDetVal=NULL)
Definition: mpr_base.cc:3060
int Typ()
Definition: subexpr.cc:976
const char * Name()
Definition: subexpr.h:121
Definition: mpr_base.h:98
#define nPrint(a)
only for debug, over any initalized currRing
Definition: numbers.h:46
void * data
Definition: subexpr.h:89
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:12
Definition: intvec.h:14
int j
Definition: myNF.cc:70
bool success()
Definition: mpr_numeric.h:162
void arrange()
Definition: mpr_numeric.cc:896
int i
Definition: cfEzgcd.cc:123
void PrintS(const char *s)
Definition: reporter.cc:294
void solve_all()
Definition: mpr_numeric.cc:871
#define IDELEMS(i)
Definition: simpleideals.h:24
mprState mprIdealCheck(const ideal theIdeal, const char *name, uResultant::resMatType mtype, BOOLEAN rmatrix=false)
rootContainer ** interpolateDenseSP(BOOLEAN matchUp=false, const number subDetVal=NULL)
Definition: mpr_base.cc:2922
#define nDelete(n)
Definition: numbers.h:16
leftv next
Definition: subexpr.h:87
static BOOLEAN rField_is_long_C(const ring r)
Definition: ring.h:491
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:38
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
#define nIsZero(n)
Definition: numbers.h:19
#define NULL
Definition: omList.c:10
static BOOLEAN rField_is_long_R(const ring r)
Definition: ring.h:488
void * Data()
Definition: subexpr.cc:1118
size_t gmp_output_digits
Definition: mpr_complex.cc:44
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:62
int BOOLEAN
Definition: auxiliary.h:131
lists listOfRoots(rootArranger *self, const unsigned int oprec)
Definition: ipshell.cc:4983
virtual number getSubDet()
Definition: mpr_base.h:37
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 4725 of file ipshell.cc.

4726 {
4727  int i;
4728  ideal p,w;
4729  p= (ideal)arg1->Data();
4730  w= (ideal)arg2->Data();
4731 
4732  // w[0] = f(p^0)
4733  // w[1] = f(p^1)
4734  // ...
4735  // p can be a vector of numbers (multivariate polynom)
4736  // or one number (univariate polynom)
4737  // tdg = deg(f)
4738 
4739  int n= IDELEMS( p );
4740  int m= IDELEMS( w );
4741  int tdg= (int)(long)arg3->Data();
4742 
4743  res->data= (void*)NULL;
4744 
4745  // check the input
4746  if ( tdg < 1 )
4747  {
4748  WerrorS("Last input parameter must be > 0!");
4749  return TRUE;
4750  }
4751  if ( n != rVar(currRing) )
4752  {
4753  Werror("Size of first input ideal must be equal to %d!",rVar(currRing));
4754  return TRUE;
4755  }
4756  if ( m != (int)pow((double)tdg+1,(double)n) )
4757  {
4758  Werror("Size of second input ideal must be equal to %d!",
4759  (int)pow((double)tdg+1,(double)n));
4760  return TRUE;
4761  }
4762  if ( !(rField_is_Q(currRing) /* ||
4763  rField_is_R() || rField_is_long_R() ||
4764  rField_is_long_C()*/ ) )
4765  {
4766  WerrorS("Ground field not implemented!");
4767  return TRUE;
4768  }
4769 
4770  number tmp;
4771  number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
4772  for ( i= 0; i < n; i++ )
4773  {
4774  pevpoint[i]=nInit(0);
4775  if ( (p->m)[i] )
4776  {
4777  tmp = pGetCoeff( (p->m)[i] );
4778  if ( nIsZero(tmp) || nIsOne(tmp) || nIsMOne(tmp) )
4779  {
4780  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4781  WerrorS("Elements of first input ideal must not be equal to -1, 0, 1!");
4782  return TRUE;
4783  }
4784  } else tmp= NULL;
4785  if ( !nIsZero(tmp) )
4786  {
4787  if ( !pIsConstant((p->m)[i]))
4788  {
4789  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4790  WerrorS("Elements of first input ideal must be numbers!");
4791  return TRUE;
4792  }
4793  pevpoint[i]= nCopy( tmp );
4794  }
4795  }
4796 
4797  number *wresults= (number *)omAlloc( m * sizeof( number ) );
4798  for ( i= 0; i < m; i++ )
4799  {
4800  wresults[i]= nInit(0);
4801  if ( (w->m)[i] && !nIsZero(pGetCoeff((w->m)[i])) )
4802  {
4803  if ( !pIsConstant((w->m)[i]))
4804  {
4805  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4806  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4807  WerrorS("Elements of second input ideal must be numbers!");
4808  return TRUE;
4809  }
4810  wresults[i]= nCopy(pGetCoeff((w->m)[i]));
4811  }
4812  }
4813 
4814  vandermonde vm( m, n, tdg, pevpoint, FALSE );
4815  number *ncpoly= vm.interpolateDense( wresults );
4816  // do not free ncpoly[]!!
4817  poly rpoly= vm.numvec2poly( ncpoly );
4818 
4819  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4820  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4821 
4822  res->data= (void*)rpoly;
4823  return FALSE;
4824 }
vandermonde system solver for interpolating polynomials from their values
Definition: mpr_numeric.h:28
#define FALSE
Definition: auxiliary.h:140
return P p
Definition: myNF.cc:203
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:537
#define TRUE
Definition: auxiliary.h:144
#define nIsOne(n)
Definition: numbers.h:25
void * ADDRESS
Definition: auxiliary.h:161
void WerrorS(const char *s)
Definition: feFopen.cc:24
#define nIsMOne(n)
Definition: numbers.h:26
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy ...
Definition: monomials.h:51
#define omAlloc(size)
Definition: omAllocDecl.h:210
void * data
Definition: subexpr.h:89
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:12
int m
Definition: cfEzgcd.cc:119
#define pIsConstant(p)
like above, except that Comp might be != 0
Definition: polys.h:209
int i
Definition: cfEzgcd.cc:123
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:458
#define IDELEMS(i)
Definition: simpleideals.h:24
#define nIsZero(n)
Definition: numbers.h:19
#define NULL
Definition: omList.c:10
const CanonicalForm & w
Definition: facAbsFact.cc:55
#define nCopy(n)
Definition: numbers.h:15
void * Data()
Definition: subexpr.cc:1118
polyrec * poly
Definition: hilb.h:10
#define nInit(i)
Definition: numbers.h:24
Rational pow(const Rational &a, int e)
Definition: GMPrat.cc:418
void Werror(const char *fmt,...)
Definition: reporter.cc:199