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

#define DEFAULT_DIGITS   30

Definition at line 13 of file mpr_inout.h.

#define MPR_DENSE   1

Definition at line 15 of file mpr_inout.h.

#define MPR_SPARSE   2

Definition at line 16 of file mpr_inout.h.

Function Documentation

BOOLEAN loNewtonP ( leftv  res,
leftv  arg1 
)

compute Newton Polytopes of input polynomials

Definition at line 4144 of file ipshell.cc.

4145 {
4146  res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
4147  return FALSE;
4148 }
#define FALSE
Definition: auxiliary.h:140
ideal loNewtonPolytope(const ideal id)
Definition: mpr_base.cc:3192
void * data
Definition: subexpr.h:89
void * Data()
Definition: subexpr.cc:1097
BOOLEAN loSimplex ( leftv  res,
leftv  args 
)

Implementation of the Simplex Algorithm.

For args, see class simplex.

Definition at line 4150 of file ipshell.cc.

4151 {
4152  if ( !(rField_is_long_R(currRing)) )
4153  {
4154  WerrorS("Ground field not implemented!");
4155  return TRUE;
4156  }
4157 
4158  simplex * LP;
4159  matrix m;
4160 
4161  leftv v= args;
4162  if ( v->Typ() != MATRIX_CMD ) // 1: matrix
4163  return TRUE;
4164  else
4165  m= (matrix)(v->CopyD());
4166 
4167  LP = new simplex(MATROWS(m),MATCOLS(m));
4168  LP->mapFromMatrix(m);
4169 
4170  v= v->next;
4171  if ( v->Typ() != INT_CMD ) // 2: m = number of constraints
4172  return TRUE;
4173  else
4174  LP->m= (int)(long)(v->Data());
4175 
4176  v= v->next;
4177  if ( v->Typ() != INT_CMD ) // 3: n = number of variables
4178  return TRUE;
4179  else
4180  LP->n= (int)(long)(v->Data());
4181 
4182  v= v->next;
4183  if ( v->Typ() != INT_CMD ) // 4: m1 = number of <= constraints
4184  return TRUE;
4185  else
4186  LP->m1= (int)(long)(v->Data());
4187 
4188  v= v->next;
4189  if ( v->Typ() != INT_CMD ) // 5: m2 = number of >= constraints
4190  return TRUE;
4191  else
4192  LP->m2= (int)(long)(v->Data());
4193 
4194  v= v->next;
4195  if ( v->Typ() != INT_CMD ) // 6: m3 = number of == constraints
4196  return TRUE;
4197  else
4198  LP->m3= (int)(long)(v->Data());
4199 
4200 #ifdef mprDEBUG_PROT
4201  Print("m (constraints) %d\n",LP->m);
4202  Print("n (columns) %d\n",LP->n);
4203  Print("m1 (<=) %d\n",LP->m1);
4204  Print("m2 (>=) %d\n",LP->m2);
4205  Print("m3 (==) %d\n",LP->m3);
4206 #endif
4207 
4208  LP->compute();
4209 
4210  lists lres= (lists)omAlloc( sizeof(slists) );
4211  lres->Init( 6 );
4212 
4213  lres->m[0].rtyp= MATRIX_CMD; // output matrix
4214  lres->m[0].data=(void*)LP->mapToMatrix(m);
4215 
4216  lres->m[1].rtyp= INT_CMD; // found a solution?
4217  lres->m[1].data=(void*)(long)LP->icase;
4218 
4219  lres->m[2].rtyp= INTVEC_CMD;
4220  lres->m[2].data=(void*)LP->posvToIV();
4221 
4222  lres->m[3].rtyp= INTVEC_CMD;
4223  lres->m[3].data=(void*)LP->zrovToIV();
4224 
4225  lres->m[4].rtyp= INT_CMD;
4226  lres->m[4].data=(void*)(long)LP->m;
4227 
4228  lres->m[5].rtyp= INT_CMD;
4229  lres->m[5].data=(void*)(long)LP->n;
4230 
4231  res->data= (void*)lres;
4232 
4233  return FALSE;
4234 }
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:85
Definition: lists.h:22
#define FALSE
Definition: auxiliary.h:140
Linear Programming / Linear Optimization using Simplex - Algorithm.
Definition: mpr_numeric.h:194
#define TRUE
Definition: auxiliary.h:144
intvec * zrovToIV()
void WerrorS(const char *s)
Definition: feFopen.cc:23
int Typ()
Definition: subexpr.cc:955
#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:12
intvec * posvToIV()
BOOLEAN mapFromMatrix(matrix m)
ip_smatrix * matrix
int m
Definition: cfEzgcd.cc:119
Definition: tok.h:88
leftv next
Definition: subexpr.h:87
INLINE_THIS void Init(int l=0)
Definition: lists.h:66
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:491
int rtyp
Definition: subexpr.h:92
void * Data()
Definition: subexpr.cc:1097
#define MATROWS(i)
Definition: matpol.h:27
int icase
Definition: mpr_numeric.h:201
void * CopyD(int t)
Definition: subexpr.cc:662
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 4259 of file ipshell.cc.

4260 {
4261 
4262  poly gls;
4263  gls= (poly)(arg1->Data());
4264  int howclean= (int)(long)arg3->Data();
4265 
4266  if ( !(rField_is_R(currRing) ||
4267  rField_is_Q(currRing) ||
4270  {
4271  WerrorS("Ground field not implemented!");
4272  return TRUE;
4273  }
4274 
4277  {
4278  unsigned long int ii = (unsigned long int)arg2->Data();
4279  setGMPFloatDigits( ii, ii );
4280  }
4281 
4282  if ( gls == NULL || pIsConstant( gls ) )
4283  {
4284  WerrorS("Input polynomial is constant!");
4285  return TRUE;
4286  }
4287 
4288  int ldummy;
4289  int deg= currRing->pLDeg( gls, &ldummy, currRing );
4290  // int deg= pDeg( gls );
4291  // int len= pLength( gls );
4292  int i,vpos=0;
4293  poly piter;
4294  lists elist;
4295  lists rlist;
4296 
4297  elist= (lists)omAlloc( sizeof(slists) );
4298  elist->Init( 0 );
4299 
4300  if ( rVar(currRing) > 1 )
4301  {
4302  piter= gls;
4303  for ( i= 1; i <= rVar(currRing); i++ )
4304  if ( pGetExp( piter, i ) )
4305  {
4306  vpos= i;
4307  break;
4308  }
4309  while ( piter )
4310  {
4311  for ( i= 1; i <= rVar(currRing); i++ )
4312  if ( (vpos != i) && (pGetExp( piter, i ) != 0) )
4313  {
4314  WerrorS("The input polynomial must be univariate!");
4315  return TRUE;
4316  }
4317  pIter( piter );
4318  }
4319  }
4320 
4321  rootContainer * roots= new rootContainer();
4322  number * pcoeffs= (number *)omAlloc( (deg+1) * sizeof( number ) );
4323  piter= gls;
4324  for ( i= deg; i >= 0; i-- )
4325  {
4326  //if ( piter ) Print("deg %d, pDeg(piter) %d\n",i,pTotaldegree(piter));
4327  if ( piter && pTotaldegree(piter) == i )
4328  {
4329  pcoeffs[i]= nCopy( pGetCoeff( piter ) );
4330  //nPrint( pcoeffs[i] );PrintS(" ");
4331  pIter( piter );
4332  }
4333  else
4334  {
4335  pcoeffs[i]= nInit(0);
4336  }
4337  }
4338 
4339 #ifdef mprDEBUG_PROT
4340  for (i=deg; i >= 0; i--)
4341  {
4342  nPrint( pcoeffs[i] );PrintS(" ");
4343  }
4344  PrintLn();
4345 #endif
4346 
4347  roots->fillContainer( pcoeffs, NULL, 1, deg, rootContainer::onepoly, 1 );
4348  roots->solver( howclean );
4349 
4350  int elem= roots->getAnzRoots();
4351  char *dummy;
4352  int j;
4353 
4354  rlist= (lists)omAlloc( sizeof(slists) );
4355  rlist->Init( elem );
4356 
4358  {
4359  for ( j= 0; j < elem; j++ )
4360  {
4361  rlist->m[j].rtyp=NUMBER_CMD;
4362  rlist->m[j].data=(void *)nCopy((number)(roots->getRoot(j)));
4363  //rlist->m[j].data=(void *)(number)(roots->getRoot(j));
4364  }
4365  }
4366  else
4367  {
4368  for ( j= 0; j < elem; j++ )
4369  {
4370  dummy = complexToStr( (*roots)[j], gmp_output_digits, currRing->cf );
4371  rlist->m[j].rtyp=STRING_CMD;
4372  rlist->m[j].data=(void *)dummy;
4373  }
4374  }
4375 
4376  elist->Clean();
4377  //omFreeSize( (ADDRESS) elist, sizeof(slists) );
4378 
4379  // this is (via fillContainer) the same data as in root
4380  //for ( i= deg; i >= 0; i-- ) nDelete( &pcoeffs[i] );
4381  //omFreeSize( (ADDRESS) pcoeffs, (deg+1) * sizeof( number ) );
4382 
4383  delete roots;
4384 
4385  res->rtyp= LIST_CMD;
4386  res->data= (void*)rlist;
4387 
4388  return FALSE;
4389 }
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:322
Definition: lists.h:22
#define FALSE
Definition: auxiliary.h:140
static BOOLEAN rField_is_R(const ring r)
Definition: ring.h:467
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:540
#define TRUE
Definition: auxiliary.h:144
bool solver(const int polishmode=PM_NONE)
Definition: mpr_numeric.cc:450
void WerrorS(const char *s)
Definition: feFopen.cc:23
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:12
#define pGetExp(p, i)
Exponent.
Definition: polys.h:41
int j
Definition: myNF.cc:70
static long pTotaldegree(poly p)
Definition: polys.h:253
#define pIsConstant(p)
like above, except that Comp might be != 0
Definition: polys.h:209
void fillContainer(number *_coeffs, number *_ievpoint, const int _var, const int _tdg, const rootType _rt, const int _anz)
Definition: mpr_numeric.cc:313
int i
Definition: cfEzgcd.cc:123
void PrintS(const char *s)
Definition: reporter.cc:294
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:461
gmp_complex * getRoot(const int i)
Definition: mpr_numeric.h:88
static BOOLEAN rField_is_long_C(const ring r)
Definition: ring.h:494
INLINE_THIS void Init(int l=0)
Definition: lists.h:66
#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:491
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:1097
Definition: tok.h:96
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
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 4236 of file ipshell.cc.

4237 {
4238  ideal gls = (ideal)(arg1->Data());
4239  int imtype= (int)(long)arg2->Data();
4240 
4241  uResultant::resMatType mtype= determineMType( imtype );
4242 
4243  // check input ideal ( = polynomial system )
4244  if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk )
4245  {
4246  return TRUE;
4247  }
4248 
4249  uResultant *resMat= new uResultant( gls, mtype, false );
4250  if (resMat!=NULL)
4251  {
4252  res->rtyp = MODUL_CMD;
4253  res->data= (void*)resMat->accessResMat()->getMatrix();
4254  if (!errorreported) delete resMat;
4255  }
4256  return errorreported;
4257 }
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:144
uResultant::resMatType determineMType(int imtype)
Definition: mpr_inout.cc:135
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)
Definition: mpr_inout.cc:94
short errorreported
Definition: feFopen.cc:22
#define NULL
Definition: omList.c:10
int rtyp
Definition: subexpr.h:92
void * Data()
Definition: subexpr.cc:1097
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 4492 of file ipshell.cc.

4493 {
4494  leftv v= args;
4495 
4496  ideal gls;
4497  int imtype;
4498  int howclean;
4499 
4500  // get ideal
4501  if ( v->Typ() != IDEAL_CMD )
4502  return TRUE;
4503  else gls= (ideal)(v->Data());
4504  v= v->next;
4505 
4506  // get resultant matrix type to use (0,1)
4507  if ( v->Typ() != INT_CMD )
4508  return TRUE;
4509  else imtype= (int)(long)v->Data();
4510  v= v->next;
4511 
4512  if (imtype==0)
4513  {
4514  ideal test_id=idInit(1,1);
4515  int j;
4516  for(j=IDELEMS(gls)-1;j>=0;j--)
4517  {
4518  if (gls->m[j]!=NULL)
4519  {
4520  test_id->m[0]=gls->m[j];
4521  intvec *dummy_w=id_QHomWeight(test_id, currRing);
4522  if (dummy_w!=NULL)
4523  {
4524  WerrorS("Newton polytope not of expected dimension");
4525  delete dummy_w;
4526  return TRUE;
4527  }
4528  }
4529  }
4530  }
4531 
4532  // get and set precision in digits ( > 0 )
4533  if ( v->Typ() != INT_CMD )
4534  return TRUE;
4535  else if ( !(rField_is_R(currRing) || rField_is_long_R(currRing) || \
4537  {
4538  unsigned long int ii=(unsigned long int)v->Data();
4539  setGMPFloatDigits( ii, ii );
4540  }
4541  v= v->next;
4542 
4543  // get interpolation steps (0,1,2)
4544  if ( v->Typ() != INT_CMD )
4545  return TRUE;
4546  else howclean= (int)(long)v->Data();
4547 
4548  uResultant::resMatType mtype= determineMType( imtype );
4549  int i,count;
4550  lists listofroots= NULL;
4551  number smv= NULL;
4552  BOOLEAN interpolate_det= (mtype==uResultant::denseResMat)?TRUE:FALSE;
4553 
4554  //emptylist= (lists)omAlloc( sizeof(slists) );
4555  //emptylist->Init( 0 );
4556 
4557  //res->rtyp = LIST_CMD;
4558  //res->data= (void *)emptylist;
4559 
4560  // check input ideal ( = polynomial system )
4561  if ( mprIdealCheck( gls, args->Name(), mtype ) != mprOk )
4562  {
4563  return TRUE;
4564  }
4565 
4566  uResultant * ures;
4567  rootContainer ** iproots;
4568  rootContainer ** muiproots;
4569  rootArranger * arranger;
4570 
4571  // main task 1: setup of resultant matrix
4572  ures= new uResultant( gls, mtype );
4573  if ( ures->accessResMat()->initState() != resMatrixBase::ready )
4574  {
4575  WerrorS("Error occurred during matrix setup!");
4576  return TRUE;
4577  }
4578 
4579  // if dense resultant, check if minor nonsingular
4580  if ( mtype == uResultant::denseResMat )
4581  {
4582  smv= ures->accessResMat()->getSubDet();
4583 #ifdef mprDEBUG_PROT
4584  PrintS("// Determinant of submatrix: ");nPrint(smv);PrintLn();
4585 #endif
4586  if ( nIsZero(smv) )
4587  {
4588  WerrorS("Unsuitable input ideal: Minor of resultant matrix is singular!");
4589  return TRUE;
4590  }
4591  }
4592 
4593  // main task 2: Interpolate specialized resultant polynomials
4594  if ( interpolate_det )
4595  iproots= ures->interpolateDenseSP( false, smv );
4596  else
4597  iproots= ures->specializeInU( false, smv );
4598 
4599  // main task 3: Interpolate specialized resultant polynomials
4600  if ( interpolate_det )
4601  muiproots= ures->interpolateDenseSP( true, smv );
4602  else
4603  muiproots= ures->specializeInU( true, smv );
4604 
4605 #ifdef mprDEBUG_PROT
4606  int c= iproots[0]->getAnzElems();
4607  for (i=0; i < c; i++) pWrite(iproots[i]->getPoly());
4608  c= muiproots[0]->getAnzElems();
4609  for (i=0; i < c; i++) pWrite(muiproots[i]->getPoly());
4610 #endif
4611 
4612  // main task 4: Compute roots of specialized polys and match them up
4613  arranger= new rootArranger( iproots, muiproots, howclean );
4614  arranger->solve_all();
4615 
4616  // get list of roots
4617  if ( arranger->success() )
4618  {
4619  arranger->arrange();
4620  listofroots= listOfRoots(arranger, gmp_output_digits );
4621  }
4622  else
4623  {
4624  WerrorS("Solver was unable to find any roots!");
4625  return TRUE;
4626  }
4627 
4628  // free everything
4629  count= iproots[0]->getAnzElems();
4630  for (i=0; i < count; i++) delete iproots[i];
4631  omFreeSize( (ADDRESS) iproots, count * sizeof(rootContainer*) );
4632  count= muiproots[0]->getAnzElems();
4633  for (i=0; i < count; i++) delete muiproots[i];
4634  omFreeSize( (ADDRESS) muiproots, count * sizeof(rootContainer*) );
4635 
4636  delete ures;
4637  delete arranger;
4638  nDelete( &smv );
4639 
4640  res->data= (void *)listofroots;
4641 
4642  //emptylist->Clean();
4643  // omFreeSize( (ADDRESS) emptylist, sizeof(slists) );
4644 
4645  return FALSE;
4646 }
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:322
Base class for solving 0-dim poly systems using u-resultant.
Definition: mpr_base.h:62
Definition: tok.h:85
Definition: lists.h:22
virtual IStateType initState() const
Definition: mpr_base.h:41
#define FALSE
Definition: auxiliary.h:140
static BOOLEAN rField_is_R(const ring r)
Definition: ring.h:467
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:144
uResultant::resMatType determineMType(int imtype)
Definition: mpr_inout.cc:135
void * ADDRESS
Definition: auxiliary.h:161
void pWrite(poly p)
Definition: polys.h:279
void WerrorS(const char *s)
Definition: feFopen.cc:23
int getAnzElems()
Definition: mpr_numeric.h:95
rootContainer ** specializeInU(BOOLEAN matchUp=false, const number subDetVal=NULL)
Definition: mpr_base.cc:3060
int Typ()
Definition: subexpr.cc:955
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:12
Definition: intvec.h:16
int j
Definition: myNF.cc:70
bool success()
Definition: mpr_numeric.h:162
void arrange()
Definition: mpr_numeric.cc:896
int i
Definition: cfEzgcd.cc:123
void PrintS(const char *s)
Definition: reporter.cc:294
void solve_all()
Definition: mpr_numeric.cc:871
#define IDELEMS(i)
Definition: simpleideals.h:24
mprState mprIdealCheck(const ideal theIdeal, const char *name, uResultant::resMatType mtype, BOOLEAN rmatrix=false)
Definition: mpr_inout.cc:94
rootContainer ** interpolateDenseSP(BOOLEAN matchUp=false, const number subDetVal=NULL)
Definition: mpr_base.cc:2922
#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:494
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:491
void * Data()
Definition: subexpr.cc:1097
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
int BOOLEAN
Definition: auxiliary.h:131
lists listOfRoots(rootArranger *self, const unsigned int oprec)
Definition: ipshell.cc:4649
virtual number getSubDet()
Definition: mpr_base.h:37
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 4391 of file ipshell.cc.

4392 {
4393  int i;
4394  ideal p,w;
4395  p= (ideal)arg1->Data();
4396  w= (ideal)arg2->Data();
4397 
4398  // w[0] = f(p^0)
4399  // w[1] = f(p^1)
4400  // ...
4401  // p can be a vector of numbers (multivariate polynom)
4402  // or one number (univariate polynom)
4403  // tdg = deg(f)
4404 
4405  int n= IDELEMS( p );
4406  int m= IDELEMS( w );
4407  int tdg= (int)(long)arg3->Data();
4408 
4409  res->data= (void*)NULL;
4410 
4411  // check the input
4412  if ( tdg < 1 )
4413  {
4414  WerrorS("Last input parameter must be > 0!");
4415  return TRUE;
4416  }
4417  if ( n != rVar(currRing) )
4418  {
4419  Werror("Size of first input ideal must be equal to %d!",rVar(currRing));
4420  return TRUE;
4421  }
4422  if ( m != (int)pow((double)tdg+1,(double)n) )
4423  {
4424  Werror("Size of second input ideal must be equal to %d!",
4425  (int)pow((double)tdg+1,(double)n));
4426  return TRUE;
4427  }
4428  if ( !(rField_is_Q(currRing) /* ||
4429  rField_is_R() || rField_is_long_R() ||
4430  rField_is_long_C()*/ ) )
4431  {
4432  WerrorS("Ground field not implemented!");
4433  return TRUE;
4434  }
4435 
4436  number tmp;
4437  number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
4438  for ( i= 0; i < n; i++ )
4439  {
4440  pevpoint[i]=nInit(0);
4441  if ( (p->m)[i] )
4442  {
4443  tmp = pGetCoeff( (p->m)[i] );
4444  if ( nIsZero(tmp) || nIsOne(tmp) || nIsMOne(tmp) )
4445  {
4446  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4447  WerrorS("Elements of first input ideal must not be equal to -1, 0, 1!");
4448  return TRUE;
4449  }
4450  } else tmp= NULL;
4451  if ( !nIsZero(tmp) )
4452  {
4453  if ( !pIsConstant((p->m)[i]))
4454  {
4455  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4456  WerrorS("Elements of first input ideal must be numbers!");
4457  return TRUE;
4458  }
4459  pevpoint[i]= nCopy( tmp );
4460  }
4461  }
4462 
4463  number *wresults= (number *)omAlloc( m * sizeof( number ) );
4464  for ( i= 0; i < m; i++ )
4465  {
4466  wresults[i]= nInit(0);
4467  if ( (w->m)[i] && !nIsZero(pGetCoeff((w->m)[i])) )
4468  {
4469  if ( !pIsConstant((w->m)[i]))
4470  {
4471  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4472  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4473  WerrorS("Elements of second input ideal must be numbers!");
4474  return TRUE;
4475  }
4476  wresults[i]= nCopy(pGetCoeff((w->m)[i]));
4477  }
4478  }
4479 
4480  vandermonde vm( m, n, tdg, pevpoint, FALSE );
4481  number *ncpoly= vm.interpolateDense( wresults );
4482  // do not free ncpoly[]!!
4483  poly rpoly= vm.numvec2poly( ncpoly );
4484 
4485  omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4486  omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4487 
4488  res->data= (void*)rpoly;
4489  return FALSE;
4490 }
vandermonde system solver for interpolating polynomials from their values
Definition: mpr_numeric.h:28
#define FALSE
Definition: auxiliary.h:140
return P p
Definition: myNF.cc:203
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
const CanonicalForm CFMap CFMap int &both_non_zero int n
Definition: cfEzgcd.cc:52
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:540
#define TRUE
Definition: auxiliary.h:144
#define nIsOne(n)
Definition: numbers.h:25
void * ADDRESS
Definition: auxiliary.h:161
void WerrorS(const char *s)
Definition: feFopen.cc:23
#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:12
int m
Definition: cfEzgcd.cc:119
#define pIsConstant(p)
like above, except that Comp might be != 0
Definition: polys.h:209
int i
Definition: cfEzgcd.cc:123
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:461
#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:1097
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:199