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

4484 {
4485  res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
4486  return FALSE;
4487 }
#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:1137

◆ loSimplex()

BOOLEAN loSimplex ( leftv  res,
leftv  args 
)

Implementation of the Simplex Algorithm.

For args, see class simplex.

Definition at line 4489 of file ipshell.cc.

4490 {
4491  if ( !(rField_is_long_R(currRing)) )
4492  {
4493  WerrorS("Ground field not implemented!");
4494  return TRUE;
4495  }
4496 
4497  simplex * LP;
4498  matrix m;
4499 
4500  leftv v= args;
4501  if ( v->Typ() != MATRIX_CMD ) // 1: matrix
4502  return TRUE;
4503  else
4504  m= (matrix)(v->CopyD());
4505 
4506  LP = new simplex(MATROWS(m),MATCOLS(m));
4507  LP->mapFromMatrix(m);
4508 
4509  v= v->next;
4510  if ( v->Typ() != INT_CMD ) // 2: m = number of constraints
4511  return TRUE;
4512  else
4513  LP->m= (int)(long)(v->Data());
4514 
4515  v= v->next;
4516  if ( v->Typ() != INT_CMD ) // 3: n = number of variables
4517  return TRUE;
4518  else
4519  LP->n= (int)(long)(v->Data());
4520 
4521  v= v->next;
4522  if ( v->Typ() != INT_CMD ) // 4: m1 = number of <= constraints
4523  return TRUE;
4524  else
4525  LP->m1= (int)(long)(v->Data());
4526 
4527  v= v->next;
4528  if ( v->Typ() != INT_CMD ) // 5: m2 = number of >= constraints
4529  return TRUE;
4530  else
4531  LP->m2= (int)(long)(v->Data());
4532 
4533  v= v->next;
4534  if ( v->Typ() != INT_CMD ) // 6: m3 = number of == constraints
4535  return TRUE;
4536  else
4537  LP->m3= (int)(long)(v->Data());
4538 
4539 #ifdef mprDEBUG_PROT
4540  Print("m (constraints) %d\n",LP->m);
4541  Print("n (columns) %d\n",LP->n);
4542  Print("m1 (<=) %d\n",LP->m1);
4543  Print("m2 (>=) %d\n",LP->m2);
4544  Print("m3 (==) %d\n",LP->m3);
4545 #endif
4546 
4547  LP->compute();
4548 
4549  lists lres= (lists)omAlloc( sizeof(slists) );
4550  lres->Init( 6 );
4551 
4552  lres->m[0].rtyp= MATRIX_CMD; // output matrix
4553  lres->m[0].data=(void*)LP->mapToMatrix(m);
4554 
4555  lres->m[1].rtyp= INT_CMD; // found a solution?
4556  lres->m[1].data=(void*)(long)LP->icase;
4557 
4558  lres->m[2].rtyp= INTVEC_CMD;
4559  lres->m[2].data=(void*)LP->posvToIV();
4560 
4561  lres->m[3].rtyp= INTVEC_CMD;
4562  lres->m[3].data=(void*)LP->zrovToIV();
4563 
4564  lres->m[4].rtyp= INT_CMD;
4565  lres->m[4].data=(void*)(long)LP->m;
4566 
4567  lres->m[5].rtyp= INT_CMD;
4568  lres->m[5].data=(void*)(long)LP->n;
4569 
4570  res->data= (void*)lres;
4571 
4572  return FALSE;
4573 }
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:995
#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:1137
#define MATROWS(i)
Definition: matpol.h:27
int icase
Definition: mpr_numeric.h:201
void * CopyD(int t)
Definition: subexpr.cc:707

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

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

4576 {
4577  ideal gls = (ideal)(arg1->Data());
4578  int imtype= (int)(long)arg2->Data();
4579 
4580  uResultant::resMatType mtype= determineMType( imtype );
4581 
4582  // check input ideal ( = polynomial system )
4583  if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk )
4584  {
4585  return TRUE;
4586  }
4587 
4588  uResultant *resMat= new uResultant( gls, mtype, false );
4589  if (resMat!=NULL)
4590  {
4591  res->rtyp = MODUL_CMD;
4592  res->data= (void*)resMat->accessResMat()->getMatrix();
4593  if (!errorreported) delete resMat;
4594  }
4595  return errorreported;
4596 }
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:1137

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

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

4728 {
4729  int i;
4730  ideal p,w;
4731  p= (ideal)arg1->Data();
4732  w= (ideal)arg2->Data();
4733 
4734  // w[0] = f(p^0)
4735  // w[1] = f(p^1)
4736  // ...
4737  // p can be a vector of numbers (multivariate polynom)
4738  // or one number (univariate polynom)
4739  // tdg = deg(f)
4740 
4741  int n= IDELEMS( p );
4742  int m= IDELEMS( w );
4743  int tdg= (int)(long)arg3->Data();
4744 
4745  res->data= (void*)NULL;
4746 
4747  // check the input
4748  if ( tdg < 1 )
4749  {
4750  WerrorS("Last input parameter must be > 0!");
4751  return TRUE;
4752  }
4753  if ( n != rVar(currRing) )
4754  {
4755  Werror("Size of first input ideal must be equal to %d!",rVar(currRing));
4756  return TRUE;
4757  }
4758  if ( m != (int)pow((double)tdg+1,(double)n) )
4759  {
4760  Werror("Size of second input ideal must be equal to %d!",
4761  (int)pow((double)tdg+1,(double)n));
4762  return TRUE;
4763  }
4764  if ( !(rField_is_Q(currRing) /* ||
4765  rField_is_R() || rField_is_long_R() ||
4766  rField_is_long_C()*/ ) )
4767  {
4768  WerrorS("Ground field not implemented!");
4769  return TRUE;
4770  }
4771 
4772  number tmp;
4773  number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
4774  for ( i= 0; i < n; i++ )
4775  {
4776  pevpoint[i]=nInit(0);
4777  if ( (p->m)[i] )
4778  {
4779  tmp = pGetCoeff( (p->m)[i] );
4780  if ( nIsZero(tmp) || nIsOne(tmp) || nIsMOne(tmp) )
4781  {
4782  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4783  WerrorS("Elements of first input ideal must not be equal to -1, 0, 1!");
4784  return TRUE;
4785  }
4786  } else tmp= NULL;
4787  if ( !nIsZero(tmp) )
4788  {
4789  if ( !pIsConstant((p->m)[i]))
4790  {
4791  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4792  WerrorS("Elements of first input ideal must be numbers!");
4793  return TRUE;
4794  }
4795  pevpoint[i]= nCopy( tmp );
4796  }
4797  }
4798 
4799  number *wresults= (number *)omAlloc( m * sizeof( number ) );
4800  for ( i= 0; i < m; i++ )
4801  {
4802  wresults[i]= nInit(0);
4803  if ( (w->m)[i] && !nIsZero(pGetCoeff((w->m)[i])) )
4804  {
4805  if ( !pIsConstant((w->m)[i]))
4806  {
4807  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4808  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4809  WerrorS("Elements of second input ideal must be numbers!");
4810  return TRUE;
4811  }
4812  wresults[i]= nCopy(pGetCoeff((w->m)[i]));
4813  }
4814  }
4815 
4816  vandermonde vm( m, n, tdg, pevpoint, FALSE );
4817  number *ncpoly= vm.interpolateDense( wresults );
4818  // do not free ncpoly[]!!
4819  poly rpoly= vm.numvec2poly( ncpoly );
4820 
4821  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4822  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4823 
4824  res->data= (void*)rpoly;
4825  return FALSE;
4826 }
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:1137
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