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

4428 {
4429  res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
4430  return FALSE;
4431 }
#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:1111
BOOLEAN loSimplex ( leftv  res,
leftv  args 
)

Implementation of the Simplex Algorithm.

For args, see class simplex.

Definition at line 4433 of file ipshell.cc.

4434 {
4435  if ( !(rField_is_long_R(currRing)) )
4436  {
4437  WerrorS("Ground field not implemented!");
4438  return TRUE;
4439  }
4440 
4441  simplex * LP;
4442  matrix m;
4443 
4444  leftv v= args;
4445  if ( v->Typ() != MATRIX_CMD ) // 1: matrix
4446  return TRUE;
4447  else
4448  m= (matrix)(v->CopyD());
4449 
4450  LP = new simplex(MATROWS(m),MATCOLS(m));
4451  LP->mapFromMatrix(m);
4452 
4453  v= v->next;
4454  if ( v->Typ() != INT_CMD ) // 2: m = number of constraints
4455  return TRUE;
4456  else
4457  LP->m= (int)(long)(v->Data());
4458 
4459  v= v->next;
4460  if ( v->Typ() != INT_CMD ) // 3: n = number of variables
4461  return TRUE;
4462  else
4463  LP->n= (int)(long)(v->Data());
4464 
4465  v= v->next;
4466  if ( v->Typ() != INT_CMD ) // 4: m1 = number of <= constraints
4467  return TRUE;
4468  else
4469  LP->m1= (int)(long)(v->Data());
4470 
4471  v= v->next;
4472  if ( v->Typ() != INT_CMD ) // 5: m2 = number of >= constraints
4473  return TRUE;
4474  else
4475  LP->m2= (int)(long)(v->Data());
4476 
4477  v= v->next;
4478  if ( v->Typ() != INT_CMD ) // 6: m3 = number of == constraints
4479  return TRUE;
4480  else
4481  LP->m3= (int)(long)(v->Data());
4482 
4483 #ifdef mprDEBUG_PROT
4484  Print("m (constraints) %d\n",LP->m);
4485  Print("n (columns) %d\n",LP->n);
4486  Print("m1 (<=) %d\n",LP->m1);
4487  Print("m2 (>=) %d\n",LP->m2);
4488  Print("m3 (==) %d\n",LP->m3);
4489 #endif
4490 
4491  LP->compute();
4492 
4493  lists lres= (lists)omAlloc( sizeof(slists) );
4494  lres->Init( 6 );
4495 
4496  lres->m[0].rtyp= MATRIX_CMD; // output matrix
4497  lres->m[0].data=(void*)LP->mapToMatrix(m);
4498 
4499  lres->m[1].rtyp= INT_CMD; // found a solution?
4500  lres->m[1].data=(void*)(long)LP->icase;
4501 
4502  lres->m[2].rtyp= INTVEC_CMD;
4503  lres->m[2].data=(void*)LP->posvToIV();
4504 
4505  lres->m[3].rtyp= INTVEC_CMD;
4506  lres->m[3].data=(void*)LP->zrovToIV();
4507 
4508  lres->m[4].rtyp= INT_CMD;
4509  lres->m[4].data=(void*)(long)LP->m;
4510 
4511  lres->m[5].rtyp= INT_CMD;
4512  lres->m[5].data=(void*)(long)LP->n;
4513 
4514  res->data= (void*)lres;
4515 
4516  return FALSE;
4517 }
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:85
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:969
#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
Definition: tok.h:88
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:491
int rtyp
Definition: subexpr.h:92
void * Data()
Definition: subexpr.cc:1111
#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 4542 of file ipshell.cc.

4543 {
4544 
4545  poly gls;
4546  gls= (poly)(arg1->Data());
4547  int howclean= (int)(long)arg3->Data();
4548 
4549  if ( !(rField_is_R(currRing) ||
4550  rField_is_Q(currRing) ||
4553  {
4554  WerrorS("Ground field not implemented!");
4555  return TRUE;
4556  }
4557 
4560  {
4561  unsigned long int ii = (unsigned long int)arg2->Data();
4562  setGMPFloatDigits( ii, ii );
4563  }
4564 
4565  if ( gls == NULL || pIsConstant( gls ) )
4566  {
4567  WerrorS("Input polynomial is constant!");
4568  return TRUE;
4569  }
4570 
4571  int ldummy;
4572  int deg= currRing->pLDeg( gls, &ldummy, currRing );
4573  // int deg= pDeg( gls );
4574  // int len= pLength( gls );
4575  int i,vpos=0;
4576  poly piter;
4577  lists elist;
4578  lists rlist;
4579 
4580  elist= (lists)omAlloc( sizeof(slists) );
4581  elist->Init( 0 );
4582 
4583  if ( rVar(currRing) > 1 )
4584  {
4585  piter= gls;
4586  for ( i= 1; i <= rVar(currRing); i++ )
4587  if ( pGetExp( piter, i ) )
4588  {
4589  vpos= i;
4590  break;
4591  }
4592  while ( piter )
4593  {
4594  for ( i= 1; i <= rVar(currRing); i++ )
4595  if ( (vpos != i) && (pGetExp( piter, i ) != 0) )
4596  {
4597  WerrorS("The input polynomial must be univariate!");
4598  return TRUE;
4599  }
4600  pIter( piter );
4601  }
4602  }
4603 
4604  rootContainer * roots= new rootContainer();
4605  number * pcoeffs= (number *)omAlloc( (deg+1) * sizeof( number ) );
4606  piter= gls;
4607  for ( i= deg; i >= 0; i-- )
4608  {
4609  //if ( piter ) Print("deg %d, pDeg(piter) %d\n",i,pTotaldegree(piter));
4610  if ( piter && pTotaldegree(piter) == i )
4611  {
4612  pcoeffs[i]= nCopy( pGetCoeff( piter ) );
4613  //nPrint( pcoeffs[i] );PrintS(" ");
4614  pIter( piter );
4615  }
4616  else
4617  {
4618  pcoeffs[i]= nInit(0);
4619  }
4620  }
4621 
4622 #ifdef mprDEBUG_PROT
4623  for (i=deg; i >= 0; i--)
4624  {
4625  nPrint( pcoeffs[i] );PrintS(" ");
4626  }
4627  PrintLn();
4628 #endif
4629 
4630  roots->fillContainer( pcoeffs, NULL, 1, deg, rootContainer::onepoly, 1 );
4631  roots->solver( howclean );
4632 
4633  int elem= roots->getAnzRoots();
4634  char *dummy;
4635  int j;
4636 
4637  rlist= (lists)omAlloc( sizeof(slists) );
4638  rlist->Init( elem );
4639 
4641  {
4642  for ( j= 0; j < elem; j++ )
4643  {
4644  rlist->m[j].rtyp=NUMBER_CMD;
4645  rlist->m[j].data=(void *)nCopy((number)(roots->getRoot(j)));
4646  //rlist->m[j].data=(void *)(number)(roots->getRoot(j));
4647  }
4648  }
4649  else
4650  {
4651  for ( j= 0; j < elem; j++ )
4652  {
4653  dummy = complexToStr( (*roots)[j], gmp_output_digits, currRing->cf );
4654  rlist->m[j].rtyp=STRING_CMD;
4655  rlist->m[j].data=(void *)dummy;
4656  }
4657  }
4658 
4659  elist->Clean();
4660  //omFreeSize( (ADDRESS) elist, sizeof(slists) );
4661 
4662  // this is (via fillContainer) the same data as in root
4663  //for ( i= deg; i >= 0; i-- ) nDelete( &pcoeffs[i] );
4664  //omFreeSize( (ADDRESS) pcoeffs, (deg+1) * sizeof( number ) );
4665 
4666  delete roots;
4667 
4668  res->rtyp= LIST_CMD;
4669  res->data= (void*)rlist;
4670 
4671  return FALSE;
4672 }
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:467
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:540
#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:461
gmp_complex * getRoot(const int i)
Definition: mpr_numeric.h:88
static BOOLEAN rField_is_long_C(const ring r)
Definition: ring.h:494
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:491
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:1111
Definition: tok.h:96
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 4519 of file ipshell.cc.

4520 {
4521  ideal gls = (ideal)(arg1->Data());
4522  int imtype= (int)(long)arg2->Data();
4523 
4524  uResultant::resMatType mtype= determineMType( imtype );
4525 
4526  // check input ideal ( = polynomial system )
4527  if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk )
4528  {
4529  return TRUE;
4530  }
4531 
4532  uResultant *resMat= new uResultant( gls, mtype, false );
4533  if (resMat!=NULL)
4534  {
4535  res->rtyp = MODUL_CMD;
4536  res->data= (void*)resMat->accessResMat()->getMatrix();
4537  if (!errorreported) delete resMat;
4538  }
4539  return errorreported;
4540 }
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:1111
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 4775 of file ipshell.cc.

4776 {
4777  leftv v= args;
4778 
4779  ideal gls;
4780  int imtype;
4781  int howclean;
4782 
4783  // get ideal
4784  if ( v->Typ() != IDEAL_CMD )
4785  return TRUE;
4786  else gls= (ideal)(v->Data());
4787  v= v->next;
4788 
4789  // get resultant matrix type to use (0,1)
4790  if ( v->Typ() != INT_CMD )
4791  return TRUE;
4792  else imtype= (int)(long)v->Data();
4793  v= v->next;
4794 
4795  if (imtype==0)
4796  {
4797  ideal test_id=idInit(1,1);
4798  int j;
4799  for(j=IDELEMS(gls)-1;j>=0;j--)
4800  {
4801  if (gls->m[j]!=NULL)
4802  {
4803  test_id->m[0]=gls->m[j];
4804  intvec *dummy_w=id_QHomWeight(test_id, currRing);
4805  if (dummy_w!=NULL)
4806  {
4807  WerrorS("Newton polytope not of expected dimension");
4808  delete dummy_w;
4809  return TRUE;
4810  }
4811  }
4812  }
4813  }
4814 
4815  // get and set precision in digits ( > 0 )
4816  if ( v->Typ() != INT_CMD )
4817  return TRUE;
4818  else if ( !(rField_is_R(currRing) || rField_is_long_R(currRing) || \
4820  {
4821  unsigned long int ii=(unsigned long int)v->Data();
4822  setGMPFloatDigits( ii, ii );
4823  }
4824  v= v->next;
4825 
4826  // get interpolation steps (0,1,2)
4827  if ( v->Typ() != INT_CMD )
4828  return TRUE;
4829  else howclean= (int)(long)v->Data();
4830 
4831  uResultant::resMatType mtype= determineMType( imtype );
4832  int i,count;
4833  lists listofroots= NULL;
4834  number smv= NULL;
4835  BOOLEAN interpolate_det= (mtype==uResultant::denseResMat)?TRUE:FALSE;
4836 
4837  //emptylist= (lists)omAlloc( sizeof(slists) );
4838  //emptylist->Init( 0 );
4839 
4840  //res->rtyp = LIST_CMD;
4841  //res->data= (void *)emptylist;
4842 
4843  // check input ideal ( = polynomial system )
4844  if ( mprIdealCheck( gls, args->Name(), mtype ) != mprOk )
4845  {
4846  return TRUE;
4847  }
4848 
4849  uResultant * ures;
4850  rootContainer ** iproots;
4851  rootContainer ** muiproots;
4852  rootArranger * arranger;
4853 
4854  // main task 1: setup of resultant matrix
4855  ures= new uResultant( gls, mtype );
4856  if ( ures->accessResMat()->initState() != resMatrixBase::ready )
4857  {
4858  WerrorS("Error occurred during matrix setup!");
4859  return TRUE;
4860  }
4861 
4862  // if dense resultant, check if minor nonsingular
4863  if ( mtype == uResultant::denseResMat )
4864  {
4865  smv= ures->accessResMat()->getSubDet();
4866 #ifdef mprDEBUG_PROT
4867  PrintS("// Determinant of submatrix: ");nPrint(smv);PrintLn();
4868 #endif
4869  if ( nIsZero(smv) )
4870  {
4871  WerrorS("Unsuitable input ideal: Minor of resultant matrix is singular!");
4872  return TRUE;
4873  }
4874  }
4875 
4876  // main task 2: Interpolate specialized resultant polynomials
4877  if ( interpolate_det )
4878  iproots= ures->interpolateDenseSP( false, smv );
4879  else
4880  iproots= ures->specializeInU( false, smv );
4881 
4882  // main task 3: Interpolate specialized resultant polynomials
4883  if ( interpolate_det )
4884  muiproots= ures->interpolateDenseSP( true, smv );
4885  else
4886  muiproots= ures->specializeInU( true, smv );
4887 
4888 #ifdef mprDEBUG_PROT
4889  int c= iproots[0]->getAnzElems();
4890  for (i=0; i < c; i++) pWrite(iproots[i]->getPoly());
4891  c= muiproots[0]->getAnzElems();
4892  for (i=0; i < c; i++) pWrite(muiproots[i]->getPoly());
4893 #endif
4894 
4895  // main task 4: Compute roots of specialized polys and match them up
4896  arranger= new rootArranger( iproots, muiproots, howclean );
4897  arranger->solve_all();
4898 
4899  // get list of roots
4900  if ( arranger->success() )
4901  {
4902  arranger->arrange();
4903  listofroots= listOfRoots(arranger, gmp_output_digits );
4904  }
4905  else
4906  {
4907  WerrorS("Solver was unable to find any roots!");
4908  return TRUE;
4909  }
4910 
4911  // free everything
4912  count= iproots[0]->getAnzElems();
4913  for (i=0; i < count; i++) delete iproots[i];
4914  omFreeSize( (ADDRESS) iproots, count * sizeof(rootContainer*) );
4915  count= muiproots[0]->getAnzElems();
4916  for (i=0; i < count; i++) delete muiproots[i];
4917  omFreeSize( (ADDRESS) muiproots, count * sizeof(rootContainer*) );
4918 
4919  delete ures;
4920  delete arranger;
4921  nDelete( &smv );
4922 
4923  res->data= (void *)listofroots;
4924 
4925  //emptylist->Clean();
4926  // omFreeSize( (ADDRESS) emptylist, sizeof(slists) );
4927 
4928  return FALSE;
4929 }
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:85
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:467
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:969
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:16
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:494
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:491
void * Data()
Definition: subexpr.cc:1111
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:4932
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 4674 of file ipshell.cc.

4675 {
4676  int i;
4677  ideal p,w;
4678  p= (ideal)arg1->Data();
4679  w= (ideal)arg2->Data();
4680 
4681  // w[0] = f(p^0)
4682  // w[1] = f(p^1)
4683  // ...
4684  // p can be a vector of numbers (multivariate polynom)
4685  // or one number (univariate polynom)
4686  // tdg = deg(f)
4687 
4688  int n= IDELEMS( p );
4689  int m= IDELEMS( w );
4690  int tdg= (int)(long)arg3->Data();
4691 
4692  res->data= (void*)NULL;
4693 
4694  // check the input
4695  if ( tdg < 1 )
4696  {
4697  WerrorS("Last input parameter must be > 0!");
4698  return TRUE;
4699  }
4700  if ( n != rVar(currRing) )
4701  {
4702  Werror("Size of first input ideal must be equal to %d!",rVar(currRing));
4703  return TRUE;
4704  }
4705  if ( m != (int)pow((double)tdg+1,(double)n) )
4706  {
4707  Werror("Size of second input ideal must be equal to %d!",
4708  (int)pow((double)tdg+1,(double)n));
4709  return TRUE;
4710  }
4711  if ( !(rField_is_Q(currRing) /* ||
4712  rField_is_R() || rField_is_long_R() ||
4713  rField_is_long_C()*/ ) )
4714  {
4715  WerrorS("Ground field not implemented!");
4716  return TRUE;
4717  }
4718 
4719  number tmp;
4720  number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
4721  for ( i= 0; i < n; i++ )
4722  {
4723  pevpoint[i]=nInit(0);
4724  if ( (p->m)[i] )
4725  {
4726  tmp = pGetCoeff( (p->m)[i] );
4727  if ( nIsZero(tmp) || nIsOne(tmp) || nIsMOne(tmp) )
4728  {
4729  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4730  WerrorS("Elements of first input ideal must not be equal to -1, 0, 1!");
4731  return TRUE;
4732  }
4733  } else tmp= NULL;
4734  if ( !nIsZero(tmp) )
4735  {
4736  if ( !pIsConstant((p->m)[i]))
4737  {
4738  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4739  WerrorS("Elements of first input ideal must be numbers!");
4740  return TRUE;
4741  }
4742  pevpoint[i]= nCopy( tmp );
4743  }
4744  }
4745 
4746  number *wresults= (number *)omAlloc( m * sizeof( number ) );
4747  for ( i= 0; i < m; i++ )
4748  {
4749  wresults[i]= nInit(0);
4750  if ( (w->m)[i] && !nIsZero(pGetCoeff((w->m)[i])) )
4751  {
4752  if ( !pIsConstant((w->m)[i]))
4753  {
4754  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4755  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4756  WerrorS("Elements of second input ideal must be numbers!");
4757  return TRUE;
4758  }
4759  wresults[i]= nCopy(pGetCoeff((w->m)[i]));
4760  }
4761  }
4762 
4763  vandermonde vm( m, n, tdg, pevpoint, FALSE );
4764  number *ncpoly= vm.interpolateDense( wresults );
4765  // do not free ncpoly[]!!
4766  poly rpoly= vm.numvec2poly( ncpoly );
4767 
4768  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4769  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4770 
4771  res->data= (void*)rpoly;
4772  return FALSE;
4773 }
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:540
#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:461
#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:1111
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