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

4454 {
4455  res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
4456  return FALSE;
4457 }
#define FALSE
Definition: auxiliary.h:94
ideal loNewtonPolytope(const ideal id)
Definition: mpr_base.cc:3190
void * data
Definition: subexpr.h:89
void * Data()
Definition: subexpr.cc:1138

§ loSimplex()

BOOLEAN loSimplex ( leftv  res,
leftv  args 
)

Implementation of the Simplex Algorithm.

For args, see class simplex.

Definition at line 4459 of file ipshell.cc.

4460 {
4461  if ( !(rField_is_long_R(currRing)) )
4462  {
4463  WerrorS("Ground field not implemented!");
4464  return TRUE;
4465  }
4466 
4467  simplex * LP;
4468  matrix m;
4469 
4470  leftv v= args;
4471  if ( v->Typ() != MATRIX_CMD ) // 1: matrix
4472  return TRUE;
4473  else
4474  m= (matrix)(v->CopyD());
4475 
4476  LP = new simplex(MATROWS(m),MATCOLS(m));
4477  LP->mapFromMatrix(m);
4478 
4479  v= v->next;
4480  if ( v->Typ() != INT_CMD ) // 2: m = number of constraints
4481  return TRUE;
4482  else
4483  LP->m= (int)(long)(v->Data());
4484 
4485  v= v->next;
4486  if ( v->Typ() != INT_CMD ) // 3: n = number of variables
4487  return TRUE;
4488  else
4489  LP->n= (int)(long)(v->Data());
4490 
4491  v= v->next;
4492  if ( v->Typ() != INT_CMD ) // 4: m1 = number of <= constraints
4493  return TRUE;
4494  else
4495  LP->m1= (int)(long)(v->Data());
4496 
4497  v= v->next;
4498  if ( v->Typ() != INT_CMD ) // 5: m2 = number of >= constraints
4499  return TRUE;
4500  else
4501  LP->m2= (int)(long)(v->Data());
4502 
4503  v= v->next;
4504  if ( v->Typ() != INT_CMD ) // 6: m3 = number of == constraints
4505  return TRUE;
4506  else
4507  LP->m3= (int)(long)(v->Data());
4508 
4509 #ifdef mprDEBUG_PROT
4510  Print("m (constraints) %d\n",LP->m);
4511  Print("n (columns) %d\n",LP->n);
4512  Print("m1 (<=) %d\n",LP->m1);
4513  Print("m2 (>=) %d\n",LP->m2);
4514  Print("m3 (==) %d\n",LP->m3);
4515 #endif
4516 
4517  LP->compute();
4518 
4519  lists lres= (lists)omAlloc( sizeof(slists) );
4520  lres->Init( 6 );
4521 
4522  lres->m[0].rtyp= MATRIX_CMD; // output matrix
4523  lres->m[0].data=(void*)LP->mapToMatrix(m);
4524 
4525  lres->m[1].rtyp= INT_CMD; // found a solution?
4526  lres->m[1].data=(void*)(long)LP->icase;
4527 
4528  lres->m[2].rtyp= INTVEC_CMD;
4529  lres->m[2].data=(void*)LP->posvToIV();
4530 
4531  lres->m[3].rtyp= INTVEC_CMD;
4532  lres->m[3].data=(void*)LP->zrovToIV();
4533 
4534  lres->m[4].rtyp= INT_CMD;
4535  lres->m[4].data=(void*)(long)LP->m;
4536 
4537  lres->m[5].rtyp= INT_CMD;
4538  lres->m[5].data=(void*)(long)LP->n;
4539 
4540  res->data= (void*)lres;
4541 
4542  return FALSE;
4543 }
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:94
Linear Programming / Linear Optimization using Simplex - Algorithm.
Definition: mpr_numeric.h:194
#define TRUE
Definition: auxiliary.h:98
intvec * zrovToIV()
void WerrorS(const char *s)
Definition: feFopen.cc:24
int Typ()
Definition: subexpr.cc:996
#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:534
int rtyp
Definition: subexpr.h:92
void * Data()
Definition: subexpr.cc:1138
#define MATROWS(i)
Definition: matpol.h:27
int icase
Definition: mpr_numeric.h:201
void * CopyD(int t)
Definition: subexpr.cc:708

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

4569 {
4570 
4571  poly gls;
4572  gls= (poly)(arg1->Data());
4573  int howclean= (int)(long)arg3->Data();
4574 
4575  if ( !(rField_is_R(currRing) ||
4576  rField_is_Q(currRing) ||
4579  {
4580  WerrorS("Ground field not implemented!");
4581  return TRUE;
4582  }
4583 
4586  {
4587  unsigned long int ii = (unsigned long int)arg2->Data();
4588  setGMPFloatDigits( ii, ii );
4589  }
4590 
4591  if ( gls == NULL || pIsConstant( gls ) )
4592  {
4593  WerrorS("Input polynomial is constant!");
4594  return TRUE;
4595  }
4596 
4597  int ldummy;
4598  int deg= currRing->pLDeg( gls, &ldummy, currRing );
4599  int i,vpos=0;
4600  poly piter;
4601  lists elist;
4602  lists rlist;
4603 
4604  elist= (lists)omAlloc( sizeof(slists) );
4605  elist->Init( 0 );
4606 
4607  if ( rVar(currRing) > 1 )
4608  {
4609  piter= gls;
4610  for ( i= 1; i <= rVar(currRing); i++ )
4611  if ( pGetExp( piter, i ) )
4612  {
4613  vpos= i;
4614  break;
4615  }
4616  while ( piter )
4617  {
4618  for ( i= 1; i <= rVar(currRing); i++ )
4619  if ( (vpos != i) && (pGetExp( piter, i ) != 0) )
4620  {
4621  WerrorS("The input polynomial must be univariate!");
4622  return TRUE;
4623  }
4624  pIter( piter );
4625  }
4626  }
4627 
4628  rootContainer * roots= new rootContainer();
4629  number * pcoeffs= (number *)omAlloc( (deg+1) * sizeof( number ) );
4630  piter= gls;
4631  for ( i= deg; i >= 0; i-- )
4632  {
4633  if ( piter && pTotaldegree(piter) == i )
4634  {
4635  pcoeffs[i]= nCopy( pGetCoeff( piter ) );
4636  //nPrint( pcoeffs[i] );PrintS(" ");
4637  pIter( piter );
4638  }
4639  else
4640  {
4641  pcoeffs[i]= nInit(0);
4642  }
4643  }
4644 
4645 #ifdef mprDEBUG_PROT
4646  for (i=deg; i >= 0; i--)
4647  {
4648  nPrint( pcoeffs[i] );PrintS(" ");
4649  }
4650  PrintLn();
4651 #endif
4652 
4653  roots->fillContainer( pcoeffs, NULL, 1, deg, rootContainer::onepoly, 1 );
4654  roots->solver( howclean );
4655 
4656  int elem= roots->getAnzRoots();
4657  char *dummy;
4658  int j;
4659 
4660  rlist= (lists)omAlloc( sizeof(slists) );
4661  rlist->Init( elem );
4662 
4664  {
4665  for ( j= 0; j < elem; j++ )
4666  {
4667  rlist->m[j].rtyp=NUMBER_CMD;
4668  rlist->m[j].data=(void *)nCopy((number)(roots->getRoot(j)));
4669  //rlist->m[j].data=(void *)(number)(roots->getRoot(j));
4670  }
4671  }
4672  else
4673  {
4674  for ( j= 0; j < elem; j++ )
4675  {
4676  dummy = complexToStr( (*roots)[j], gmp_output_digits, currRing->cf );
4677  rlist->m[j].rtyp=STRING_CMD;
4678  rlist->m[j].data=(void *)dummy;
4679  }
4680  }
4681 
4682  elist->Clean();
4683  //omFreeSize( (ADDRESS) elist, sizeof(slists) );
4684 
4685  // this is (via fillContainer) the same data as in root
4686  //for ( i= deg; i >= 0; i-- ) nDelete( &pcoeffs[i] );
4687  //omFreeSize( (ADDRESS) pcoeffs, (deg+1) * sizeof( number ) );
4688 
4689  delete roots;
4690 
4691  res->rtyp= LIST_CMD;
4692  res->data= (void*)rlist;
4693 
4694  return FALSE;
4695 }
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:94
static BOOLEAN rField_is_R(const ring r)
Definition: ring.h:510
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:583
#define TRUE
Definition: auxiliary.h:98
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:264
#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:537
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:534
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:1138
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 4545 of file ipshell.cc.

4546 {
4547  ideal gls = (ideal)(arg1->Data());
4548  int imtype= (int)(long)arg2->Data();
4549 
4550  uResultant::resMatType mtype= determineMType( imtype );
4551 
4552  // check input ideal ( = polynomial system )
4553  if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk )
4554  {
4555  return TRUE;
4556  }
4557 
4558  uResultant *resMat= new uResultant( gls, mtype, false );
4559  if (resMat!=NULL)
4560  {
4561  res->rtyp = MODUL_CMD;
4562  res->data= (void*)resMat->accessResMat()->getMatrix();
4563  if (!errorreported) delete resMat;
4564  }
4565  return errorreported;
4566 }
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:98
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:1138

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

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

4698 {
4699  int i;
4700  ideal p,w;
4701  p= (ideal)arg1->Data();
4702  w= (ideal)arg2->Data();
4703 
4704  // w[0] = f(p^0)
4705  // w[1] = f(p^1)
4706  // ...
4707  // p can be a vector of numbers (multivariate polynom)
4708  // or one number (univariate polynom)
4709  // tdg = deg(f)
4710 
4711  int n= IDELEMS( p );
4712  int m= IDELEMS( w );
4713  int tdg= (int)(long)arg3->Data();
4714 
4715  res->data= (void*)NULL;
4716 
4717  // check the input
4718  if ( tdg < 1 )
4719  {
4720  WerrorS("Last input parameter must be > 0!");
4721  return TRUE;
4722  }
4723  if ( n != rVar(currRing) )
4724  {
4725  Werror("Size of first input ideal must be equal to %d!",rVar(currRing));
4726  return TRUE;
4727  }
4728  if ( m != (int)pow((double)tdg+1,(double)n) )
4729  {
4730  Werror("Size of second input ideal must be equal to %d!",
4731  (int)pow((double)tdg+1,(double)n));
4732  return TRUE;
4733  }
4734  if ( !(rField_is_Q(currRing) /* ||
4735  rField_is_R() || rField_is_long_R() ||
4736  rField_is_long_C()*/ ) )
4737  {
4738  WerrorS("Ground field not implemented!");
4739  return TRUE;
4740  }
4741 
4742  number tmp;
4743  number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
4744  for ( i= 0; i < n; i++ )
4745  {
4746  pevpoint[i]=nInit(0);
4747  if ( (p->m)[i] )
4748  {
4749  tmp = pGetCoeff( (p->m)[i] );
4750  if ( nIsZero(tmp) || nIsOne(tmp) || nIsMOne(tmp) )
4751  {
4752  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4753  WerrorS("Elements of first input ideal must not be equal to -1, 0, 1!");
4754  return TRUE;
4755  }
4756  } else tmp= NULL;
4757  if ( !nIsZero(tmp) )
4758  {
4759  if ( !pIsConstant((p->m)[i]))
4760  {
4761  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4762  WerrorS("Elements of first input ideal must be numbers!");
4763  return TRUE;
4764  }
4765  pevpoint[i]= nCopy( tmp );
4766  }
4767  }
4768 
4769  number *wresults= (number *)omAlloc( m * sizeof( number ) );
4770  for ( i= 0; i < m; i++ )
4771  {
4772  wresults[i]= nInit(0);
4773  if ( (w->m)[i] && !nIsZero(pGetCoeff((w->m)[i])) )
4774  {
4775  if ( !pIsConstant((w->m)[i]))
4776  {
4777  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4778  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4779  WerrorS("Elements of second input ideal must be numbers!");
4780  return TRUE;
4781  }
4782  wresults[i]= nCopy(pGetCoeff((w->m)[i]));
4783  }
4784  }
4785 
4786  vandermonde vm( m, n, tdg, pevpoint, FALSE );
4787  number *ncpoly= vm.interpolateDense( wresults );
4788  // do not free ncpoly[]!!
4789  poly rpoly= vm.numvec2poly( ncpoly );
4790 
4791  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4792  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4793 
4794  res->data= (void*)rpoly;
4795  return FALSE;
4796 }
vandermonde system solver for interpolating polynomials from their values
Definition: mpr_numeric.h:28
#define FALSE
Definition: auxiliary.h:94
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:583
#define TRUE
Definition: auxiliary.h:98
#define nIsOne(n)
Definition: numbers.h:25
void * ADDRESS
Definition: auxiliary.h:115
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:1138
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