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

4469 {
4470  res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
4471  return FALSE;
4472 }
#define FALSE
Definition: auxiliary.h:97
ideal loNewtonPolytope(const ideal id)
Definition: mpr_base.cc:3190
void * data
Definition: subexpr.h:90
void * Data()
Definition: subexpr.cc:1121

§ loSimplex()

BOOLEAN loSimplex ( leftv  res,
leftv  args 
)

Implementation of the Simplex Algorithm.

For args, see class simplex.

Definition at line 4474 of file ipshell.cc.

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

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

4584 {
4585 
4586  poly gls;
4587  gls= (poly)(arg1->Data());
4588  int howclean= (int)(long)arg3->Data();
4589 
4590  if ( !(rField_is_R(currRing) ||
4591  rField_is_Q(currRing) ||
4594  {
4595  WerrorS("Ground field not implemented!");
4596  return TRUE;
4597  }
4598 
4601  {
4602  unsigned long int ii = (unsigned long int)arg2->Data();
4603  setGMPFloatDigits( ii, ii );
4604  }
4605 
4606  if ( gls == NULL || pIsConstant( gls ) )
4607  {
4608  WerrorS("Input polynomial is constant!");
4609  return TRUE;
4610  }
4611 
4612  int ldummy;
4613  int deg= currRing->pLDeg( gls, &ldummy, currRing );
4614  int i,vpos=0;
4615  poly piter;
4616  lists elist;
4617  lists rlist;
4618 
4619  elist= (lists)omAlloc( sizeof(slists) );
4620  elist->Init( 0 );
4621 
4622  if ( rVar(currRing) > 1 )
4623  {
4624  piter= gls;
4625  for ( i= 1; i <= rVar(currRing); i++ )
4626  if ( pGetExp( piter, i ) )
4627  {
4628  vpos= i;
4629  break;
4630  }
4631  while ( piter )
4632  {
4633  for ( i= 1; i <= rVar(currRing); i++ )
4634  if ( (vpos != i) && (pGetExp( piter, i ) != 0) )
4635  {
4636  WerrorS("The input polynomial must be univariate!");
4637  return TRUE;
4638  }
4639  pIter( piter );
4640  }
4641  }
4642 
4643  rootContainer * roots= new rootContainer();
4644  number * pcoeffs= (number *)omAlloc( (deg+1) * sizeof( number ) );
4645  piter= gls;
4646  for ( i= deg; i >= 0; i-- )
4647  {
4648  if ( piter && pTotaldegree(piter) == i )
4649  {
4650  pcoeffs[i]= nCopy( pGetCoeff( piter ) );
4651  //nPrint( pcoeffs[i] );PrintS(" ");
4652  pIter( piter );
4653  }
4654  else
4655  {
4656  pcoeffs[i]= nInit(0);
4657  }
4658  }
4659 
4660 #ifdef mprDEBUG_PROT
4661  for (i=deg; i >= 0; i--)
4662  {
4663  nPrint( pcoeffs[i] );PrintS(" ");
4664  }
4665  PrintLn();
4666 #endif
4667 
4668  roots->fillContainer( pcoeffs, NULL, 1, deg, rootContainer::onepoly, 1 );
4669  roots->solver( howclean );
4670 
4671  int elem= roots->getAnzRoots();
4672  char *dummy;
4673  int j;
4674 
4675  rlist= (lists)omAlloc( sizeof(slists) );
4676  rlist->Init( elem );
4677 
4679  {
4680  for ( j= 0; j < elem; j++ )
4681  {
4682  rlist->m[j].rtyp=NUMBER_CMD;
4683  rlist->m[j].data=(void *)nCopy((number)(roots->getRoot(j)));
4684  //rlist->m[j].data=(void *)(number)(roots->getRoot(j));
4685  }
4686  }
4687  else
4688  {
4689  for ( j= 0; j < elem; j++ )
4690  {
4691  dummy = complexToStr( (*roots)[j], gmp_output_digits, currRing->cf );
4692  rlist->m[j].rtyp=STRING_CMD;
4693  rlist->m[j].data=(void *)dummy;
4694  }
4695  }
4696 
4697  elist->Clean();
4698  //omFreeSize( (ADDRESS) elist, sizeof(slists) );
4699 
4700  // this is (via fillContainer) the same data as in root
4701  //for ( i= deg; i >= 0; i-- ) nDelete( &pcoeffs[i] );
4702  //omFreeSize( (ADDRESS) pcoeffs, (deg+1) * sizeof( number ) );
4703 
4704  delete roots;
4705 
4706  res->rtyp= LIST_CMD;
4707  res->data= (void*)rlist;
4708 
4709  return FALSE;
4710 }
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:97
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:101
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:90
#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:93
#define nCopy(n)
Definition: numbers.h:15
void Clean(ring r=currRing)
Definition: lists.h:25
void * Data()
Definition: subexpr.cc:1121
Definition: tok.h:116
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

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

4561 {
4562  ideal gls = (ideal)(arg1->Data());
4563  int imtype= (int)(long)arg2->Data();
4564 
4565  uResultant::resMatType mtype= determineMType( imtype );
4566 
4567  // check input ideal ( = polynomial system )
4568  if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk )
4569  {
4570  return TRUE;
4571  }
4572 
4573  uResultant *resMat= new uResultant( gls, mtype, false );
4574  if (resMat!=NULL)
4575  {
4576  res->rtyp = MODUL_CMD;
4577  res->data= (void*)resMat->accessResMat()->getMatrix();
4578  if (!errorreported) delete resMat;
4579  }
4580  return errorreported;
4581 }
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:101
uResultant::resMatType determineMType(int imtype)
const char * Name()
Definition: subexpr.h:122
Definition: mpr_base.h:98
void * data
Definition: subexpr.h:90
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:93
void * Data()
Definition: subexpr.cc:1121

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

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

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