My Project  UNKNOWN_GIT_VERSION
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 4532 of file ipshell.cc.

4533 {
4534  res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
4535  return FALSE;
4536 }
#define FALSE
Definition: auxiliary.h:94
ideal loNewtonPolytope(const ideal id)
Definition: mpr_base.cc:3192
CanonicalForm res
Definition: facAbsFact.cc:64
void * Data()
Definition: subexpr.cc:1182

◆ loSimplex()

BOOLEAN loSimplex ( leftv  res,
leftv  args 
)

Implementation of the Simplex Algorithm.

For args, see class simplex.

Definition at line 4538 of file ipshell.cc.

4539 {
4540  if ( !(rField_is_long_R(currRing)) )
4541  {
4542  WerrorS("Ground field not implemented!");
4543  return TRUE;
4544  }
4545 
4546  simplex * LP;
4547  matrix m;
4548 
4549  leftv v= args;
4550  if ( v->Typ() != MATRIX_CMD ) // 1: matrix
4551  return TRUE;
4552  else
4553  m= (matrix)(v->CopyD());
4554 
4555  LP = new simplex(MATROWS(m),MATCOLS(m));
4556  LP->mapFromMatrix(m);
4557 
4558  v= v->next;
4559  if ( v->Typ() != INT_CMD ) // 2: m = number of constraints
4560  return TRUE;
4561  else
4562  LP->m= (int)(long)(v->Data());
4563 
4564  v= v->next;
4565  if ( v->Typ() != INT_CMD ) // 3: n = number of variables
4566  return TRUE;
4567  else
4568  LP->n= (int)(long)(v->Data());
4569 
4570  v= v->next;
4571  if ( v->Typ() != INT_CMD ) // 4: m1 = number of <= constraints
4572  return TRUE;
4573  else
4574  LP->m1= (int)(long)(v->Data());
4575 
4576  v= v->next;
4577  if ( v->Typ() != INT_CMD ) // 5: m2 = number of >= constraints
4578  return TRUE;
4579  else
4580  LP->m2= (int)(long)(v->Data());
4581 
4582  v= v->next;
4583  if ( v->Typ() != INT_CMD ) // 6: m3 = number of == constraints
4584  return TRUE;
4585  else
4586  LP->m3= (int)(long)(v->Data());
4587 
4588 #ifdef mprDEBUG_PROT
4589  Print("m (constraints) %d\n",LP->m);
4590  Print("n (columns) %d\n",LP->n);
4591  Print("m1 (<=) %d\n",LP->m1);
4592  Print("m2 (>=) %d\n",LP->m2);
4593  Print("m3 (==) %d\n",LP->m3);
4594 #endif
4595 
4596  LP->compute();
4597 
4598  lists lres= (lists)omAlloc( sizeof(slists) );
4599  lres->Init( 6 );
4600 
4601  lres->m[0].rtyp= MATRIX_CMD; // output matrix
4602  lres->m[0].data=(void*)LP->mapToMatrix(m);
4603 
4604  lres->m[1].rtyp= INT_CMD; // found a solution?
4605  lres->m[1].data=(void*)(long)LP->icase;
4606 
4607  lres->m[2].rtyp= INTVEC_CMD;
4608  lres->m[2].data=(void*)LP->posvToIV();
4609 
4610  lres->m[3].rtyp= INTVEC_CMD;
4611  lres->m[3].data=(void*)LP->zrovToIV();
4612 
4613  lres->m[4].rtyp= INT_CMD;
4614  lres->m[4].data=(void*)(long)LP->m;
4615 
4616  lres->m[5].rtyp= INT_CMD;
4617  lres->m[5].data=(void*)(long)LP->n;
4618 
4619  res->data= (void*)lres;
4620 
4621  return FALSE;
4622 }
sleftv * m
Definition: lists.h:45
Class used for (list of) interpreter objects.
Definition: subexpr.h:82
matrix mapToMatrix(matrix m)
void compute()
#define Print
Definition: emacs.cc:80
Definition: tok.h:96
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
#define omAlloc(size)
Definition: omAllocDecl.h:210
void * data
Definition: subexpr.h:88
intvec * posvToIV()
CanonicalForm res
Definition: facAbsFact.cc:64
BOOLEAN mapFromMatrix(matrix m)
int m
Definition: cfEzgcd.cc:121
Variable next() const
Definition: factory.h:137
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:27
slists * lists
Definition: mpr_numeric.h:146
static BOOLEAN rField_is_long_R(const ring r)
Definition: ring.h:533
int rtyp
Definition: subexpr.h:91
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
#define MATROWS(i)
Definition: matpol.h:26
int icase
Definition: mpr_numeric.h:201
ip_smatrix * matrix
Definition: matpol.h:31

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

4648 {
4649 
4650  poly gls;
4651  gls= (poly)(arg1->Data());
4652  int howclean= (int)(long)arg3->Data();
4653 
4654  if ( !(rField_is_R(currRing) ||
4655  rField_is_Q(currRing) ||
4658  {
4659  WerrorS("Ground field not implemented!");
4660  return TRUE;
4661  }
4662 
4665  {
4666  unsigned long int ii = (unsigned long int)arg2->Data();
4667  setGMPFloatDigits( ii, ii );
4668  }
4669 
4670  if ( gls == NULL || pIsConstant( gls ) )
4671  {
4672  WerrorS("Input polynomial is constant!");
4673  return TRUE;
4674  }
4675 
4676  int ldummy;
4677  int deg= currRing->pLDeg( gls, &ldummy, currRing );
4678  int i,vpos=0;
4679  poly piter;
4680  lists elist;
4681  lists rlist;
4682 
4683  elist= (lists)omAlloc( sizeof(slists) );
4684  elist->Init( 0 );
4685 
4686  if ( rVar(currRing) > 1 )
4687  {
4688  piter= gls;
4689  for ( i= 1; i <= rVar(currRing); i++ )
4690  if ( pGetExp( piter, i ) )
4691  {
4692  vpos= i;
4693  break;
4694  }
4695  while ( piter )
4696  {
4697  for ( i= 1; i <= rVar(currRing); i++ )
4698  if ( (vpos != i) && (pGetExp( piter, i ) != 0) )
4699  {
4700  WerrorS("The input polynomial must be univariate!");
4701  return TRUE;
4702  }
4703  pIter( piter );
4704  }
4705  }
4706 
4707  rootContainer * roots= new rootContainer();
4708  number * pcoeffs= (number *)omAlloc( (deg+1) * sizeof( number ) );
4709  piter= gls;
4710  for ( i= deg; i >= 0; i-- )
4711  {
4712  if ( piter && pTotaldegree(piter) == i )
4713  {
4714  pcoeffs[i]= nCopy( pGetCoeff( piter ) );
4715  //nPrint( pcoeffs[i] );PrintS(" ");
4716  pIter( piter );
4717  }
4718  else
4719  {
4720  pcoeffs[i]= nInit(0);
4721  }
4722  }
4723 
4724 #ifdef mprDEBUG_PROT
4725  for (i=deg; i >= 0; i--)
4726  {
4727  nPrint( pcoeffs[i] );PrintS(" ");
4728  }
4729  PrintLn();
4730 #endif
4731 
4732  roots->fillContainer( pcoeffs, NULL, 1, deg, rootContainer::onepoly, 1 );
4733  roots->solver( howclean );
4734 
4735  int elem= roots->getAnzRoots();
4736  char *dummy;
4737  int j;
4738 
4739  rlist= (lists)omAlloc( sizeof(slists) );
4740  rlist->Init( elem );
4741 
4743  {
4744  for ( j= 0; j < elem; j++ )
4745  {
4746  rlist->m[j].rtyp=NUMBER_CMD;
4747  rlist->m[j].data=(void *)nCopy((number)(roots->getRoot(j)));
4748  //rlist->m[j].data=(void *)(number)(roots->getRoot(j));
4749  }
4750  }
4751  else
4752  {
4753  for ( j= 0; j < elem; j++ )
4754  {
4755  dummy = complexToStr( (*roots)[j], gmp_output_digits, currRing->cf );
4756  rlist->m[j].rtyp=STRING_CMD;
4757  rlist->m[j].data=(void *)dummy;
4758  }
4759  }
4760 
4761  elist->Clean();
4762  //omFreeSize( (ADDRESS) elist, sizeof(slists) );
4763 
4764  // this is (via fillContainer) the same data as in root
4765  //for ( i= deg; i >= 0; i-- ) nDelete( &pcoeffs[i] );
4766  //omFreeSize( (ADDRESS) pcoeffs, (deg+1) * sizeof( number ) );
4767 
4768  delete roots;
4769 
4770  res->rtyp= LIST_CMD;
4771  res->data= (void*)rlist;
4772 
4773  return FALSE;
4774 }
complex root finder for univariate polynomials based on laguers algorithm
Definition: mpr_numeric.h:65
sleftv * m
Definition: lists.h:45
int j
Definition: facHensel.cc:105
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:509
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:582
#define TRUE
Definition: auxiliary.h:98
bool solver(const int polishmode=PM_NONE)
Definition: mpr_numeric.cc:441
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:45
#define omAlloc(size)
Definition: omAllocDecl.h:210
#define nPrint(a)
only for debug, over any initalized currRing
Definition: numbers.h:47
void * data
Definition: subexpr.h:88
#define pIter(p)
Definition: monomials.h:38
#define pGetExp(p, i)
Exponent.
Definition: polys.h:41
CanonicalForm res
Definition: facAbsFact.cc:64
static long pTotaldegree(poly p)
Definition: polys.h:276
#define pIsConstant(p)
like above, except that Comp might be != 0
Definition: polys.h:233
void fillContainer(number *_coeffs, number *_ievpoint, const int _var, const int _tdg, const rootType _rt, const int _anz)
Definition: mpr_numeric.cc:304
int i
Definition: cfEzgcd.cc:125
void PrintS(const char *s)
Definition: reporter.cc:284
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:497
gmp_complex * getRoot(const int i)
Definition: mpr_numeric.h:88
static BOOLEAN rField_is_long_C(const ring r)
Definition: ring.h:536
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:533
int rtyp
Definition: subexpr.h:91
#define nCopy(n)
Definition: numbers.h:16
void Clean(ring r=currRing)
Definition: lists.h:25
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
void * Data()
Definition: subexpr.cc:1182
Definition: tok.h:118
char * complexToStr(gmp_complex &c, const unsigned int oprec, const coeffs src)
Definition: mpr_complex.cc:705
size_t gmp_output_digits
Definition: mpr_complex.cc:43
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:61
#define nInit(i)
Definition: numbers.h:25

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

4625 {
4626  ideal gls = (ideal)(arg1->Data());
4627  int imtype= (int)(long)arg2->Data();
4628 
4629  uResultant::resMatType mtype= determineMType( imtype );
4630 
4631  // check input ideal ( = polynomial system )
4632  if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk )
4633  {
4634  return TRUE;
4635  }
4636 
4637  uResultant *resMat= new uResultant( gls, mtype, false );
4638  if (resMat!=NULL)
4639  {
4640  res->rtyp = MODUL_CMD;
4641  res->data= (void*)resMat->accessResMat()->getMatrix();
4642  if (!errorreported) delete resMat;
4643  }
4644  return errorreported;
4645 }
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:120
Definition: mpr_base.h:98
CanonicalForm res
Definition: facAbsFact.cc:64
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
void * Data()
Definition: subexpr.cc:1182

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

4878 {
4879  leftv v= args;
4880 
4881  ideal gls;
4882  int imtype;
4883  int howclean;
4884 
4885  // get ideal
4886  if ( v->Typ() != IDEAL_CMD )
4887  return TRUE;
4888  else gls= (ideal)(v->Data());
4889  v= v->next;
4890 
4891  // get resultant matrix type to use (0,1)
4892  if ( v->Typ() != INT_CMD )
4893  return TRUE;
4894  else imtype= (int)(long)v->Data();
4895  v= v->next;
4896 
4897  if (imtype==0)
4898  {
4899  ideal test_id=idInit(1,1);
4900  int j;
4901  for(j=IDELEMS(gls)-1;j>=0;j--)
4902  {
4903  if (gls->m[j]!=NULL)
4904  {
4905  test_id->m[0]=gls->m[j];
4906  intvec *dummy_w=id_QHomWeight(test_id, currRing);
4907  if (dummy_w!=NULL)
4908  {
4909  WerrorS("Newton polytope not of expected dimension");
4910  delete dummy_w;
4911  return TRUE;
4912  }
4913  }
4914  }
4915  }
4916 
4917  // get and set precision in digits ( > 0 )
4918  if ( v->Typ() != INT_CMD )
4919  return TRUE;
4920  else if ( !(rField_is_R(currRing) || rField_is_long_R(currRing) || \
4922  {
4923  unsigned long int ii=(unsigned long int)v->Data();
4924  setGMPFloatDigits( ii, ii );
4925  }
4926  v= v->next;
4927 
4928  // get interpolation steps (0,1,2)
4929  if ( v->Typ() != INT_CMD )
4930  return TRUE;
4931  else howclean= (int)(long)v->Data();
4932 
4933  uResultant::resMatType mtype= determineMType( imtype );
4934  int i,count;
4935  lists listofroots= NULL;
4936  number smv= NULL;
4937  BOOLEAN interpolate_det= (mtype==uResultant::denseResMat)?TRUE:FALSE;
4938 
4939  //emptylist= (lists)omAlloc( sizeof(slists) );
4940  //emptylist->Init( 0 );
4941 
4942  //res->rtyp = LIST_CMD;
4943  //res->data= (void *)emptylist;
4944 
4945  // check input ideal ( = polynomial system )
4946  if ( mprIdealCheck( gls, args->Name(), mtype ) != mprOk )
4947  {
4948  return TRUE;
4949  }
4950 
4951  uResultant * ures;
4952  rootContainer ** iproots;
4953  rootContainer ** muiproots;
4954  rootArranger * arranger;
4955 
4956  // main task 1: setup of resultant matrix
4957  ures= new uResultant( gls, mtype );
4958  if ( ures->accessResMat()->initState() != resMatrixBase::ready )
4959  {
4960  WerrorS("Error occurred during matrix setup!");
4961  return TRUE;
4962  }
4963 
4964  // if dense resultant, check if minor nonsingular
4965  if ( mtype == uResultant::denseResMat )
4966  {
4967  smv= ures->accessResMat()->getSubDet();
4968 #ifdef mprDEBUG_PROT
4969  PrintS("// Determinant of submatrix: ");nPrint(smv);PrintLn();
4970 #endif
4971  if ( nIsZero(smv) )
4972  {
4973  WerrorS("Unsuitable input ideal: Minor of resultant matrix is singular!");
4974  return TRUE;
4975  }
4976  }
4977 
4978  // main task 2: Interpolate specialized resultant polynomials
4979  if ( interpolate_det )
4980  iproots= ures->interpolateDenseSP( false, smv );
4981  else
4982  iproots= ures->specializeInU( false, smv );
4983 
4984  // main task 3: Interpolate specialized resultant polynomials
4985  if ( interpolate_det )
4986  muiproots= ures->interpolateDenseSP( true, smv );
4987  else
4988  muiproots= ures->specializeInU( true, smv );
4989 
4990 #ifdef mprDEBUG_PROT
4991  int c= iproots[0]->getAnzElems();
4992  for (i=0; i < c; i++) pWrite(iproots[i]->getPoly());
4993  c= muiproots[0]->getAnzElems();
4994  for (i=0; i < c; i++) pWrite(muiproots[i]->getPoly());
4995 #endif
4996 
4997  // main task 4: Compute roots of specialized polys and match them up
4998  arranger= new rootArranger( iproots, muiproots, howclean );
4999  arranger->solve_all();
5000 
5001  // get list of roots
5002  if ( arranger->success() )
5003  {
5004  arranger->arrange();
5005  listofroots= listOfRoots(arranger, gmp_output_digits );
5006  }
5007  else
5008  {
5009  WerrorS("Solver was unable to find any roots!");
5010  return TRUE;
5011  }
5012 
5013  // free everything
5014  count= iproots[0]->getAnzElems();
5015  for (i=0; i < count; i++) delete iproots[i];
5016  omFreeSize( (ADDRESS) iproots, count * sizeof(rootContainer*) );
5017  count= muiproots[0]->getAnzElems();
5018  for (i=0; i < count; i++) delete muiproots[i];
5019  omFreeSize( (ADDRESS) muiproots, count * sizeof(rootContainer*) );
5020 
5021  delete ures;
5022  delete arranger;
5023  nDelete( &smv );
5024 
5025  res->data= (void *)listofroots;
5026 
5027  //emptylist->Clean();
5028  // omFreeSize( (ADDRESS) emptylist, sizeof(slists) );
5029 
5030  return FALSE;
5031 }
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
int j
Definition: facHensel.cc:105
Class used for (list of) interpreter objects.
Definition: subexpr.h:82
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:96
Definition: lists.h:22
#define FALSE
Definition: auxiliary.h:94
static BOOLEAN rField_is_R(const ring r)
Definition: ring.h:509
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:133
void pWrite(poly p)
Definition: polys.h:302
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:3061
const char * Name()
Definition: subexpr.h:120
Definition: mpr_base.h:98
#define nPrint(a)
only for debug, over any initalized currRing
Definition: numbers.h:47
Definition: intvec.h:17
CanonicalForm res
Definition: facAbsFact.cc:64
bool success()
Definition: mpr_numeric.h:162
void arrange()
Definition: mpr_numeric.cc:887
int i
Definition: cfEzgcd.cc:125
void PrintS(const char *s)
Definition: reporter.cc:284
void solve_all()
Definition: mpr_numeric.cc:862
Variable next() const
Definition: factory.h:137
#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:2923
#define nDelete(n)
Definition: numbers.h:17
static BOOLEAN rField_is_long_C(const ring r)
Definition: ring.h:536
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:37
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
#define nIsZero(n)
Definition: numbers.h:20
#define NULL
Definition: omList.c:10
static BOOLEAN rField_is_long_R(const ring r)
Definition: ring.h:533
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
size_t gmp_output_digits
Definition: mpr_complex.cc:43
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:61
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:5034
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 4776 of file ipshell.cc.

4777 {
4778  int i;
4779  ideal p,w;
4780  p= (ideal)arg1->Data();
4781  w= (ideal)arg2->Data();
4782 
4783  // w[0] = f(p^0)
4784  // w[1] = f(p^1)
4785  // ...
4786  // p can be a vector of numbers (multivariate polynom)
4787  // or one number (univariate polynom)
4788  // tdg = deg(f)
4789 
4790  int n= IDELEMS( p );
4791  int m= IDELEMS( w );
4792  int tdg= (int)(long)arg3->Data();
4793 
4794  res->data= (void*)NULL;
4795 
4796  // check the input
4797  if ( tdg < 1 )
4798  {
4799  WerrorS("Last input parameter must be > 0!");
4800  return TRUE;
4801  }
4802  if ( n != rVar(currRing) )
4803  {
4804  Werror("Size of first input ideal must be equal to %d!",rVar(currRing));
4805  return TRUE;
4806  }
4807  if ( m != (int)pow((double)tdg+1,(double)n) )
4808  {
4809  Werror("Size of second input ideal must be equal to %d!",
4810  (int)pow((double)tdg+1,(double)n));
4811  return TRUE;
4812  }
4813  if ( !(rField_is_Q(currRing) /* ||
4814  rField_is_R() || rField_is_long_R() ||
4815  rField_is_long_C()*/ ) )
4816  {
4817  WerrorS("Ground field not implemented!");
4818  return TRUE;
4819  }
4820 
4821  number tmp;
4822  number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
4823  for ( i= 0; i < n; i++ )
4824  {
4825  pevpoint[i]=nInit(0);
4826  if ( (p->m)[i] )
4827  {
4828  tmp = pGetCoeff( (p->m)[i] );
4829  if ( nIsZero(tmp) || nIsOne(tmp) || nIsMOne(tmp) )
4830  {
4831  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4832  WerrorS("Elements of first input ideal must not be equal to -1, 0, 1!");
4833  return TRUE;
4834  }
4835  } else tmp= NULL;
4836  if ( !nIsZero(tmp) )
4837  {
4838  if ( !pIsConstant((p->m)[i]))
4839  {
4840  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4841  WerrorS("Elements of first input ideal must be numbers!");
4842  return TRUE;
4843  }
4844  pevpoint[i]= nCopy( tmp );
4845  }
4846  }
4847 
4848  number *wresults= (number *)omAlloc( m * sizeof( number ) );
4849  for ( i= 0; i < m; i++ )
4850  {
4851  wresults[i]= nInit(0);
4852  if ( (w->m)[i] && !nIsZero(pGetCoeff((w->m)[i])) )
4853  {
4854  if ( !pIsConstant((w->m)[i]))
4855  {
4856  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4857  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4858  WerrorS("Elements of second input ideal must be numbers!");
4859  return TRUE;
4860  }
4861  wresults[i]= nCopy(pGetCoeff((w->m)[i]));
4862  }
4863  }
4864 
4865  vandermonde vm( m, n, tdg, pevpoint, FALSE );
4866  number *ncpoly= vm.interpolateDense( wresults );
4867  // do not free ncpoly[]!!
4868  poly rpoly= vm.numvec2poly( ncpoly );
4869 
4870  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4871  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4872 
4873  res->data= (void*)rpoly;
4874  return FALSE;
4875 }
vandermonde system solver for interpolating polynomials from their values
Definition: mpr_numeric.h:28
#define FALSE
Definition: auxiliary.h:94
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:582
#define TRUE
Definition: auxiliary.h:98
#define nIsOne(n)
Definition: numbers.h:26
void * ADDRESS
Definition: auxiliary.h:133
void WerrorS(const char *s)
Definition: feFopen.cc:24
#define nIsMOne(n)
Definition: numbers.h:27
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:45
#define omAlloc(size)
Definition: omAllocDecl.h:210
CanonicalForm res
Definition: facAbsFact.cc:64
int m
Definition: cfEzgcd.cc:121
#define pIsConstant(p)
like above, except that Comp might be != 0
Definition: polys.h:233
int i
Definition: cfEzgcd.cc:125
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:497
#define IDELEMS(i)
Definition: simpleideals.h:24
#define nIsZero(n)
Definition: numbers.h:20
#define NULL
Definition: omList.c:10
const CanonicalForm & w
Definition: facAbsFact.cc:55
#define nCopy(n)
Definition: numbers.h:16
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
void * Data()
Definition: subexpr.cc:1182
int p
Definition: cfModGcd.cc:4019
#define nInit(i)
Definition: numbers.h:25
Rational pow(const Rational &a, int e)
Definition: GMPrat.cc:414
void Werror(const char *fmt,...)
Definition: reporter.cc:189