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

4473 {
4474  res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
4475  return FALSE;
4476 }
#define FALSE
Definition: auxiliary.h:95
ideal loNewtonPolytope(const ideal id)
Definition: mpr_base.cc:3190
void * data
Definition: subexpr.h:89
void * Data()
Definition: subexpr.cc:1146

§ loSimplex()

BOOLEAN loSimplex ( leftv  res,
leftv  args 
)

Implementation of the Simplex Algorithm.

For args, see class simplex.

Definition at line 4478 of file ipshell.cc.

4479 {
4480  if ( !(rField_is_long_R(currRing)) )
4481  {
4482  WerrorS("Ground field not implemented!");
4483  return TRUE;
4484  }
4485 
4486  simplex * LP;
4487  matrix m;
4488 
4489  leftv v= args;
4490  if ( v->Typ() != MATRIX_CMD ) // 1: matrix
4491  return TRUE;
4492  else
4493  m= (matrix)(v->CopyD());
4494 
4495  LP = new simplex(MATROWS(m),MATCOLS(m));
4496  LP->mapFromMatrix(m);
4497 
4498  v= v->next;
4499  if ( v->Typ() != INT_CMD ) // 2: m = number of constraints
4500  return TRUE;
4501  else
4502  LP->m= (int)(long)(v->Data());
4503 
4504  v= v->next;
4505  if ( v->Typ() != INT_CMD ) // 3: n = number of variables
4506  return TRUE;
4507  else
4508  LP->n= (int)(long)(v->Data());
4509 
4510  v= v->next;
4511  if ( v->Typ() != INT_CMD ) // 4: m1 = number of <= constraints
4512  return TRUE;
4513  else
4514  LP->m1= (int)(long)(v->Data());
4515 
4516  v= v->next;
4517  if ( v->Typ() != INT_CMD ) // 5: m2 = number of >= constraints
4518  return TRUE;
4519  else
4520  LP->m2= (int)(long)(v->Data());
4521 
4522  v= v->next;
4523  if ( v->Typ() != INT_CMD ) // 6: m3 = number of == constraints
4524  return TRUE;
4525  else
4526  LP->m3= (int)(long)(v->Data());
4527 
4528 #ifdef mprDEBUG_PROT
4529  Print("m (constraints) %d\n",LP->m);
4530  Print("n (columns) %d\n",LP->n);
4531  Print("m1 (<=) %d\n",LP->m1);
4532  Print("m2 (>=) %d\n",LP->m2);
4533  Print("m3 (==) %d\n",LP->m3);
4534 #endif
4535 
4536  LP->compute();
4537 
4538  lists lres= (lists)omAlloc( sizeof(slists) );
4539  lres->Init( 6 );
4540 
4541  lres->m[0].rtyp= MATRIX_CMD; // output matrix
4542  lres->m[0].data=(void*)LP->mapToMatrix(m);
4543 
4544  lres->m[1].rtyp= INT_CMD; // found a solution?
4545  lres->m[1].data=(void*)(long)LP->icase;
4546 
4547  lres->m[2].rtyp= INTVEC_CMD;
4548  lres->m[2].data=(void*)LP->posvToIV();
4549 
4550  lres->m[3].rtyp= INTVEC_CMD;
4551  lres->m[3].data=(void*)LP->zrovToIV();
4552 
4553  lres->m[4].rtyp= INT_CMD;
4554  lres->m[4].data=(void*)(long)LP->m;
4555 
4556  lres->m[5].rtyp= INT_CMD;
4557  lres->m[5].data=(void*)(long)LP->n;
4558 
4559  res->data= (void*)lres;
4560 
4561  return FALSE;
4562 }
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:95
Definition: lists.h:22
#define FALSE
Definition: auxiliary.h:95
Linear Programming / Linear Optimization using Simplex - Algorithm.
Definition: mpr_numeric.h:194
#define TRUE
Definition: auxiliary.h:99
intvec * zrovToIV()
void WerrorS(const char *s)
Definition: feFopen.cc:24
int Typ()
Definition: subexpr.cc:1004
#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:10
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)
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:531
int rtyp
Definition: subexpr.h:92
void * Data()
Definition: subexpr.cc:1146
#define MATROWS(i)
Definition: matpol.h:27
int icase
Definition: mpr_numeric.h:201
void * CopyD(int t)
Definition: subexpr.cc:714

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

4588 {
4589 
4590  poly gls;
4591  gls= (poly)(arg1->Data());
4592  int howclean= (int)(long)arg3->Data();
4593 
4594  if ( !(rField_is_R(currRing) ||
4595  rField_is_Q(currRing) ||
4598  {
4599  WerrorS("Ground field not implemented!");
4600  return TRUE;
4601  }
4602 
4605  {
4606  unsigned long int ii = (unsigned long int)arg2->Data();
4607  setGMPFloatDigits( ii, ii );
4608  }
4609 
4610  if ( gls == NULL || pIsConstant( gls ) )
4611  {
4612  WerrorS("Input polynomial is constant!");
4613  return TRUE;
4614  }
4615 
4616  int ldummy;
4617  int deg= currRing->pLDeg( gls, &ldummy, currRing );
4618  int i,vpos=0;
4619  poly piter;
4620  lists elist;
4621  lists rlist;
4622 
4623  elist= (lists)omAlloc( sizeof(slists) );
4624  elist->Init( 0 );
4625 
4626  if ( rVar(currRing) > 1 )
4627  {
4628  piter= gls;
4629  for ( i= 1; i <= rVar(currRing); i++ )
4630  if ( pGetExp( piter, i ) )
4631  {
4632  vpos= i;
4633  break;
4634  }
4635  while ( piter )
4636  {
4637  for ( i= 1; i <= rVar(currRing); i++ )
4638  if ( (vpos != i) && (pGetExp( piter, i ) != 0) )
4639  {
4640  WerrorS("The input polynomial must be univariate!");
4641  return TRUE;
4642  }
4643  pIter( piter );
4644  }
4645  }
4646 
4647  rootContainer * roots= new rootContainer();
4648  number * pcoeffs= (number *)omAlloc( (deg+1) * sizeof( number ) );
4649  piter= gls;
4650  for ( i= deg; i >= 0; i-- )
4651  {
4652  if ( piter && pTotaldegree(piter) == i )
4653  {
4654  pcoeffs[i]= nCopy( pGetCoeff( piter ) );
4655  //nPrint( pcoeffs[i] );PrintS(" ");
4656  pIter( piter );
4657  }
4658  else
4659  {
4660  pcoeffs[i]= nInit(0);
4661  }
4662  }
4663 
4664 #ifdef mprDEBUG_PROT
4665  for (i=deg; i >= 0; i--)
4666  {
4667  nPrint( pcoeffs[i] );PrintS(" ");
4668  }
4669  PrintLn();
4670 #endif
4671 
4672  roots->fillContainer( pcoeffs, NULL, 1, deg, rootContainer::onepoly, 1 );
4673  roots->solver( howclean );
4674 
4675  int elem= roots->getAnzRoots();
4676  char *dummy;
4677  int j;
4678 
4679  rlist= (lists)omAlloc( sizeof(slists) );
4680  rlist->Init( elem );
4681 
4683  {
4684  for ( j= 0; j < elem; j++ )
4685  {
4686  rlist->m[j].rtyp=NUMBER_CMD;
4687  rlist->m[j].data=(void *)nCopy((number)(roots->getRoot(j)));
4688  //rlist->m[j].data=(void *)(number)(roots->getRoot(j));
4689  }
4690  }
4691  else
4692  {
4693  for ( j= 0; j < elem; j++ )
4694  {
4695  dummy = complexToStr( (*roots)[j], gmp_output_digits, currRing->cf );
4696  rlist->m[j].rtyp=STRING_CMD;
4697  rlist->m[j].data=(void *)dummy;
4698  }
4699  }
4700 
4701  elist->Clean();
4702  //omFreeSize( (ADDRESS) elist, sizeof(slists) );
4703 
4704  // this is (via fillContainer) the same data as in root
4705  //for ( i= deg; i >= 0; i-- ) nDelete( &pcoeffs[i] );
4706  //omFreeSize( (ADDRESS) pcoeffs, (deg+1) * sizeof( number ) );
4707 
4708  delete roots;
4709 
4710  res->rtyp= LIST_CMD;
4711  res->data= (void*)rlist;
4712 
4713  return FALSE;
4714 }
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:310
Definition: lists.h:22
#define FALSE
Definition: auxiliary.h:95
static BOOLEAN rField_is_R(const ring r)
Definition: ring.h:507
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:580
#define TRUE
Definition: auxiliary.h:99
bool solver(const int polishmode=PM_NONE)
Definition: mpr_numeric.cc:449
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:10
#define pGetExp(p, i)
Exponent.
Definition: polys.h:41
int j
Definition: myNF.cc:70
static long pTotaldegree(poly p)
Definition: polys.h:265
#define pIsConstant(p)
like above, except that Comp might be != 0
Definition: polys.h:221
void fillContainer(number *_coeffs, number *_ievpoint, const int _var, const int _tdg, const rootType _rt, const int _anz)
Definition: mpr_numeric.cc:312
int i
Definition: cfEzgcd.cc:123
void PrintS(const char *s)
Definition: reporter.cc:284
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:501
gmp_complex * getRoot(const int i)
Definition: mpr_numeric.h:88
static BOOLEAN rField_is_long_C(const ring r)
Definition: ring.h:534
INLINE_THIS void Init(int l=0)
#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:531
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:1146
Definition: tok.h:117
char * complexToStr(gmp_complex &c, const unsigned int oprec, const coeffs src)
Definition: mpr_complex.cc:706
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

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

4565 {
4566  ideal gls = (ideal)(arg1->Data());
4567  int imtype= (int)(long)arg2->Data();
4568 
4569  uResultant::resMatType mtype= determineMType( imtype );
4570 
4571  // check input ideal ( = polynomial system )
4572  if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk )
4573  {
4574  return TRUE;
4575  }
4576 
4577  uResultant *resMat= new uResultant( gls, mtype, false );
4578  if (resMat!=NULL)
4579  {
4580  res->rtyp = MODUL_CMD;
4581  res->data= (void*)resMat->accessResMat()->getMatrix();
4582  if (!errorreported) delete resMat;
4583  }
4584  return errorreported;
4585 }
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:99
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:1146

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

4818 {
4819  leftv v= args;
4820 
4821  ideal gls;
4822  int imtype;
4823  int howclean;
4824 
4825  // get ideal
4826  if ( v->Typ() != IDEAL_CMD )
4827  return TRUE;
4828  else gls= (ideal)(v->Data());
4829  v= v->next;
4830 
4831  // get resultant matrix type to use (0,1)
4832  if ( v->Typ() != INT_CMD )
4833  return TRUE;
4834  else imtype= (int)(long)v->Data();
4835  v= v->next;
4836 
4837  if (imtype==0)
4838  {
4839  ideal test_id=idInit(1,1);
4840  int j;
4841  for(j=IDELEMS(gls)-1;j>=0;j--)
4842  {
4843  if (gls->m[j]!=NULL)
4844  {
4845  test_id->m[0]=gls->m[j];
4846  intvec *dummy_w=id_QHomWeight(test_id, currRing);
4847  if (dummy_w!=NULL)
4848  {
4849  WerrorS("Newton polytope not of expected dimension");
4850  delete dummy_w;
4851  return TRUE;
4852  }
4853  }
4854  }
4855  }
4856 
4857  // get and set precision in digits ( > 0 )
4858  if ( v->Typ() != INT_CMD )
4859  return TRUE;
4860  else if ( !(rField_is_R(currRing) || rField_is_long_R(currRing) || \
4862  {
4863  unsigned long int ii=(unsigned long int)v->Data();
4864  setGMPFloatDigits( ii, ii );
4865  }
4866  v= v->next;
4867 
4868  // get interpolation steps (0,1,2)
4869  if ( v->Typ() != INT_CMD )
4870  return TRUE;
4871  else howclean= (int)(long)v->Data();
4872 
4873  uResultant::resMatType mtype= determineMType( imtype );
4874  int i,count;
4875  lists listofroots= NULL;
4876  number smv= NULL;
4877  BOOLEAN interpolate_det= (mtype==uResultant::denseResMat)?TRUE:FALSE;
4878 
4879  //emptylist= (lists)omAlloc( sizeof(slists) );
4880  //emptylist->Init( 0 );
4881 
4882  //res->rtyp = LIST_CMD;
4883  //res->data= (void *)emptylist;
4884 
4885  // check input ideal ( = polynomial system )
4886  if ( mprIdealCheck( gls, args->Name(), mtype ) != mprOk )
4887  {
4888  return TRUE;
4889  }
4890 
4891  uResultant * ures;
4892  rootContainer ** iproots;
4893  rootContainer ** muiproots;
4894  rootArranger * arranger;
4895 
4896  // main task 1: setup of resultant matrix
4897  ures= new uResultant( gls, mtype );
4898  if ( ures->accessResMat()->initState() != resMatrixBase::ready )
4899  {
4900  WerrorS("Error occurred during matrix setup!");
4901  return TRUE;
4902  }
4903 
4904  // if dense resultant, check if minor nonsingular
4905  if ( mtype == uResultant::denseResMat )
4906  {
4907  smv= ures->accessResMat()->getSubDet();
4908 #ifdef mprDEBUG_PROT
4909  PrintS("// Determinant of submatrix: ");nPrint(smv);PrintLn();
4910 #endif
4911  if ( nIsZero(smv) )
4912  {
4913  WerrorS("Unsuitable input ideal: Minor of resultant matrix is singular!");
4914  return TRUE;
4915  }
4916  }
4917 
4918  // main task 2: Interpolate specialized resultant polynomials
4919  if ( interpolate_det )
4920  iproots= ures->interpolateDenseSP( false, smv );
4921  else
4922  iproots= ures->specializeInU( false, smv );
4923 
4924  // main task 3: Interpolate specialized resultant polynomials
4925  if ( interpolate_det )
4926  muiproots= ures->interpolateDenseSP( true, smv );
4927  else
4928  muiproots= ures->specializeInU( true, smv );
4929 
4930 #ifdef mprDEBUG_PROT
4931  int c= iproots[0]->getAnzElems();
4932  for (i=0; i < c; i++) pWrite(iproots[i]->getPoly());
4933  c= muiproots[0]->getAnzElems();
4934  for (i=0; i < c; i++) pWrite(muiproots[i]->getPoly());
4935 #endif
4936 
4937  // main task 4: Compute roots of specialized polys and match them up
4938  arranger= new rootArranger( iproots, muiproots, howclean );
4939  arranger->solve_all();
4940 
4941  // get list of roots
4942  if ( arranger->success() )
4943  {
4944  arranger->arrange();
4945  listofroots= listOfRoots(arranger, gmp_output_digits );
4946  }
4947  else
4948  {
4949  WerrorS("Solver was unable to find any roots!");
4950  return TRUE;
4951  }
4952 
4953  // free everything
4954  count= iproots[0]->getAnzElems();
4955  for (i=0; i < count; i++) delete iproots[i];
4956  omFreeSize( (ADDRESS) iproots, count * sizeof(rootContainer*) );
4957  count= muiproots[0]->getAnzElems();
4958  for (i=0; i < count; i++) delete muiproots[i];
4959  omFreeSize( (ADDRESS) muiproots, count * sizeof(rootContainer*) );
4960 
4961  delete ures;
4962  delete arranger;
4963  nDelete( &smv );
4964 
4965  res->data= (void *)listofroots;
4966 
4967  //emptylist->Clean();
4968  // omFreeSize( (ADDRESS) emptylist, sizeof(slists) );
4969 
4970  return FALSE;
4971 }
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:310
Base class for solving 0-dim poly systems using u-resultant.
Definition: mpr_base.h:62
Definition: tok.h:95
Definition: lists.h:22
#define FALSE
Definition: auxiliary.h:95
static BOOLEAN rField_is_R(const ring r)
Definition: ring.h:507
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:99
uResultant::resMatType determineMType(int imtype)
void * ADDRESS
Definition: auxiliary.h:116
void pWrite(poly p)
Definition: polys.h:291
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:3059
int Typ()
Definition: subexpr.cc:1004
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:10
Definition: intvec.h:14
int j
Definition: myNF.cc:70
bool success()
Definition: mpr_numeric.h:162
void arrange()
Definition: mpr_numeric.cc:895
int i
Definition: cfEzgcd.cc:123
void PrintS(const char *s)
Definition: reporter.cc:284
void solve_all()
Definition: mpr_numeric.cc:870
#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:2921
#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:534
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:531
void * Data()
Definition: subexpr.cc:1146
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
virtual IStateType initState() const
Definition: mpr_base.h:41
int BOOLEAN
Definition: auxiliary.h:86
lists listOfRoots(rootArranger *self, const unsigned int oprec)
Definition: ipshell.cc:4974
virtual number getSubDet()
Definition: mpr_base.h:37

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

4717 {
4718  int i;
4719  ideal p,w;
4720  p= (ideal)arg1->Data();
4721  w= (ideal)arg2->Data();
4722 
4723  // w[0] = f(p^0)
4724  // w[1] = f(p^1)
4725  // ...
4726  // p can be a vector of numbers (multivariate polynom)
4727  // or one number (univariate polynom)
4728  // tdg = deg(f)
4729 
4730  int n= IDELEMS( p );
4731  int m= IDELEMS( w );
4732  int tdg= (int)(long)arg3->Data();
4733 
4734  res->data= (void*)NULL;
4735 
4736  // check the input
4737  if ( tdg < 1 )
4738  {
4739  WerrorS("Last input parameter must be > 0!");
4740  return TRUE;
4741  }
4742  if ( n != rVar(currRing) )
4743  {
4744  Werror("Size of first input ideal must be equal to %d!",rVar(currRing));
4745  return TRUE;
4746  }
4747  if ( m != (int)pow((double)tdg+1,(double)n) )
4748  {
4749  Werror("Size of second input ideal must be equal to %d!",
4750  (int)pow((double)tdg+1,(double)n));
4751  return TRUE;
4752  }
4753  if ( !(rField_is_Q(currRing) /* ||
4754  rField_is_R() || rField_is_long_R() ||
4755  rField_is_long_C()*/ ) )
4756  {
4757  WerrorS("Ground field not implemented!");
4758  return TRUE;
4759  }
4760 
4761  number tmp;
4762  number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
4763  for ( i= 0; i < n; i++ )
4764  {
4765  pevpoint[i]=nInit(0);
4766  if ( (p->m)[i] )
4767  {
4768  tmp = pGetCoeff( (p->m)[i] );
4769  if ( nIsZero(tmp) || nIsOne(tmp) || nIsMOne(tmp) )
4770  {
4771  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4772  WerrorS("Elements of first input ideal must not be equal to -1, 0, 1!");
4773  return TRUE;
4774  }
4775  } else tmp= NULL;
4776  if ( !nIsZero(tmp) )
4777  {
4778  if ( !pIsConstant((p->m)[i]))
4779  {
4780  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4781  WerrorS("Elements of first input ideal must be numbers!");
4782  return TRUE;
4783  }
4784  pevpoint[i]= nCopy( tmp );
4785  }
4786  }
4787 
4788  number *wresults= (number *)omAlloc( m * sizeof( number ) );
4789  for ( i= 0; i < m; i++ )
4790  {
4791  wresults[i]= nInit(0);
4792  if ( (w->m)[i] && !nIsZero(pGetCoeff((w->m)[i])) )
4793  {
4794  if ( !pIsConstant((w->m)[i]))
4795  {
4796  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4797  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4798  WerrorS("Elements of second input ideal must be numbers!");
4799  return TRUE;
4800  }
4801  wresults[i]= nCopy(pGetCoeff((w->m)[i]));
4802  }
4803  }
4804 
4805  vandermonde vm( m, n, tdg, pevpoint, FALSE );
4806  number *ncpoly= vm.interpolateDense( wresults );
4807  // do not free ncpoly[]!!
4808  poly rpoly= vm.numvec2poly( ncpoly );
4809 
4810  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4811  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4812 
4813  res->data= (void*)rpoly;
4814  return FALSE;
4815 }
vandermonde system solver for interpolating polynomials from their values
Definition: mpr_numeric.h:28
#define FALSE
Definition: auxiliary.h:95
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:580
#define TRUE
Definition: auxiliary.h:99
#define nIsOne(n)
Definition: numbers.h:25
void * ADDRESS
Definition: auxiliary.h:116
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:10
int m
Definition: cfEzgcd.cc:119
#define pIsConstant(p)
like above, except that Comp might be != 0
Definition: polys.h:221
int i
Definition: cfEzgcd.cc:123
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:501
#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:1146
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:189