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

4563{
4564 res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
4565 return FALSE;
4566}
#define FALSE
Definition: auxiliary.h:96
void * Data()
Definition: subexpr.cc:1154
CanonicalForm res
Definition: facAbsFact.cc:60
ideal loNewtonPolytope(const ideal id)
Definition: mpr_base.cc:3190

◆ loSimplex()

BOOLEAN loSimplex ( leftv  res,
leftv  args 
)

Implementation of the Simplex Algorithm.

For args, see class simplex.

Definition at line 4568 of file ipshell.cc.

4569{
4570 if ( !(rField_is_long_R(currRing)) )
4571 {
4572 WerrorS("Ground field not implemented!");
4573 return TRUE;
4574 }
4575
4576 simplex * LP;
4577 matrix m;
4578
4579 leftv v= args;
4580 if ( v->Typ() != MATRIX_CMD ) // 1: matrix
4581 return TRUE;
4582 else
4583 m= (matrix)(v->CopyD());
4584
4585 LP = new simplex(MATROWS(m),MATCOLS(m));
4586 LP->mapFromMatrix(m);
4587
4588 v= v->next;
4589 if ( v->Typ() != INT_CMD ) // 2: m = number of constraints
4590 return TRUE;
4591 else
4592 LP->m= (int)(long)(v->Data());
4593
4594 v= v->next;
4595 if ( v->Typ() != INT_CMD ) // 3: n = number of variables
4596 return TRUE;
4597 else
4598 LP->n= (int)(long)(v->Data());
4599
4600 v= v->next;
4601 if ( v->Typ() != INT_CMD ) // 4: m1 = number of <= constraints
4602 return TRUE;
4603 else
4604 LP->m1= (int)(long)(v->Data());
4605
4606 v= v->next;
4607 if ( v->Typ() != INT_CMD ) // 5: m2 = number of >= constraints
4608 return TRUE;
4609 else
4610 LP->m2= (int)(long)(v->Data());
4611
4612 v= v->next;
4613 if ( v->Typ() != INT_CMD ) // 6: m3 = number of == constraints
4614 return TRUE;
4615 else
4616 LP->m3= (int)(long)(v->Data());
4617
4618#ifdef mprDEBUG_PROT
4619 Print("m (constraints) %d\n",LP->m);
4620 Print("n (columns) %d\n",LP->n);
4621 Print("m1 (<=) %d\n",LP->m1);
4622 Print("m2 (>=) %d\n",LP->m2);
4623 Print("m3 (==) %d\n",LP->m3);
4624#endif
4625
4626 LP->compute();
4627
4628 lists lres= (lists)omAlloc( sizeof(slists) );
4629 lres->Init( 6 );
4630
4631 lres->m[0].rtyp= MATRIX_CMD; // output matrix
4632 lres->m[0].data=(void*)LP->mapToMatrix(m);
4633
4634 lres->m[1].rtyp= INT_CMD; // found a solution?
4635 lres->m[1].data=(void*)(long)LP->icase;
4636
4637 lres->m[2].rtyp= INTVEC_CMD;
4638 lres->m[2].data=(void*)LP->posvToIV();
4639
4640 lres->m[3].rtyp= INTVEC_CMD;
4641 lres->m[3].data=(void*)LP->zrovToIV();
4642
4643 lres->m[4].rtyp= INT_CMD;
4644 lres->m[4].data=(void*)(long)LP->m;
4645
4646 lres->m[5].rtyp= INT_CMD;
4647 lres->m[5].data=(void*)(long)LP->n;
4648
4649 res->data= (void*)lres;
4650
4651 return FALSE;
4652}
#define TRUE
Definition: auxiliary.h:100
int m
Definition: cfEzgcd.cc:128
Variable next() const
Definition: factory.h:146
Linear Programming / Linear Optimization using Simplex - Algorithm.
Definition: mpr_numeric.h:195
intvec * zrovToIV()
BOOLEAN mapFromMatrix(matrix m)
int icase
Definition: mpr_numeric.h:201
void compute()
matrix mapToMatrix(matrix m)
intvec * posvToIV()
Class used for (list of) interpreter objects.
Definition: subexpr.h:83
int rtyp
Definition: subexpr.h:91
void * data
Definition: subexpr.h:88
Definition: lists.h:24
sleftv * m
Definition: lists.h:46
INLINE_THIS void Init(int l=0)
#define Print
Definition: emacs.cc:80
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
void WerrorS(const char *s)
Definition: feFopen.cc:24
@ MATRIX_CMD
Definition: grammar.cc:286
ip_smatrix * matrix
Definition: matpol.h:43
#define MATROWS(i)
Definition: matpol.h:26
#define MATCOLS(i)
Definition: matpol.h:27
slists * lists
Definition: mpr_numeric.h:146
#define omAlloc(size)
Definition: omAllocDecl.h:210
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
static BOOLEAN rField_is_long_R(const ring r)
Definition: ring.h:543
@ INTVEC_CMD
Definition: tok.h:101
@ INT_CMD
Definition: tok.h:96

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

4678{
4679 poly gls;
4680 gls= (poly)(arg1->Data());
4681 int howclean= (int)(long)arg3->Data();
4682
4683 if ( gls == NULL || pIsConstant( gls ) )
4684 {
4685 WerrorS("Input polynomial is constant!");
4686 return TRUE;
4687 }
4688
4690 {
4691 int* r=Zp_roots(gls, currRing);
4692 lists rlist;
4693 rlist= (lists)omAlloc( sizeof(slists) );
4694 rlist->Init( r[0] );
4695 for(int i=r[0];i>0;i--)
4696 {
4697 rlist->m[i-1].data=n_Init(r[i],currRing->cf);
4698 rlist->m[i-1].rtyp=NUMBER_CMD;
4699 }
4700 omFree(r);
4701 res->data=rlist;
4702 res->rtyp= LIST_CMD;
4703 return FALSE;
4704 }
4705 if ( !(rField_is_R(currRing) ||
4709 {
4710 WerrorS("Ground field not implemented!");
4711 return TRUE;
4712 }
4713
4716 {
4717 unsigned long int ii = (unsigned long int)arg2->Data();
4718 setGMPFloatDigits( ii, ii );
4719 }
4720
4721 int ldummy;
4722 int deg= currRing->pLDeg( gls, &ldummy, currRing );
4723 int i,vpos=0;
4724 poly piter;
4725 lists elist;
4726
4727 elist= (lists)omAlloc( sizeof(slists) );
4728 elist->Init( 0 );
4729
4730 if ( rVar(currRing) > 1 )
4731 {
4732 piter= gls;
4733 for ( i= 1; i <= rVar(currRing); i++ )
4734 if ( pGetExp( piter, i ) )
4735 {
4736 vpos= i;
4737 break;
4738 }
4739 while ( piter )
4740 {
4741 for ( i= 1; i <= rVar(currRing); i++ )
4742 if ( (vpos != i) && (pGetExp( piter, i ) != 0) )
4743 {
4744 WerrorS("The input polynomial must be univariate!");
4745 return TRUE;
4746 }
4747 pIter( piter );
4748 }
4749 }
4750
4751 rootContainer * roots= new rootContainer();
4752 number * pcoeffs= (number *)omAlloc( (deg+1) * sizeof( number ) );
4753 piter= gls;
4754 for ( i= deg; i >= 0; i-- )
4755 {
4756 if ( piter && pTotaldegree(piter) == i )
4757 {
4758 pcoeffs[i]= nCopy( pGetCoeff( piter ) );
4759 //nPrint( pcoeffs[i] );PrintS(" ");
4760 pIter( piter );
4761 }
4762 else
4763 {
4764 pcoeffs[i]= nInit(0);
4765 }
4766 }
4767
4768#ifdef mprDEBUG_PROT
4769 for (i=deg; i >= 0; i--)
4770 {
4771 nPrint( pcoeffs[i] );PrintS(" ");
4772 }
4773 PrintLn();
4774#endif
4775
4776 roots->fillContainer( pcoeffs, NULL, 1, deg, rootContainer::onepoly, 1 );
4777 roots->solver( howclean );
4778
4779 int elem= roots->getAnzRoots();
4780 char *dummy;
4781 int j;
4782
4783 lists rlist;
4784 rlist= (lists)omAlloc( sizeof(slists) );
4785 rlist->Init( elem );
4786
4788 {
4789 for ( j= 0; j < elem; j++ )
4790 {
4791 rlist->m[j].rtyp=NUMBER_CMD;
4792 rlist->m[j].data=(void *)nCopy((number)(roots->getRoot(j)));
4793 //rlist->m[j].data=(void *)(number)(roots->getRoot(j));
4794 }
4795 }
4796 else
4797 {
4798 for ( j= 0; j < elem; j++ )
4799 {
4800 dummy = complexToStr( (*roots)[j], gmp_output_digits, currRing->cf );
4801 rlist->m[j].rtyp=STRING_CMD;
4802 rlist->m[j].data=(void *)dummy;
4803 }
4804 }
4805
4806 elist->Clean();
4807 //omFreeSize( (ADDRESS) elist, sizeof(slists) );
4808
4809 // this is (via fillContainer) the same data as in root
4810 //for ( i= deg; i >= 0; i-- ) nDelete( &pcoeffs[i] );
4811 //omFreeSize( (ADDRESS) pcoeffs, (deg+1) * sizeof( number ) );
4812
4813 delete roots;
4814
4815 res->data= (void*)rlist;
4816
4817 return FALSE;
4818}
int i
Definition: cfEzgcd.cc:132
int * Zp_roots(poly p, const ring r)
Definition: clapsing.cc:2188
complex root finder for univariate polynomials based on laguers algorithm
Definition: mpr_numeric.h:66
gmp_complex * getRoot(const int i)
Definition: mpr_numeric.h:88
void fillContainer(number *_coeffs, number *_ievpoint, const int _var, const int _tdg, const rootType _rt, const int _anz)
Definition: mpr_numeric.cc:300
int getAnzRoots()
Definition: mpr_numeric.h:97
bool solver(const int polishmode=PM_NONE)
Definition: mpr_numeric.cc:437
void Clean(ring r=currRing)
Definition: lists.h:26
static FORCE_INLINE number n_Init(long i, const coeffs r)
a number representing i in the given coeff field/ring r
Definition: coeffs.h:538
int j
Definition: facHensel.cc:110
@ NUMBER_CMD
Definition: grammar.cc:288
#define pIter(p)
Definition: monomials.h:37
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy
Definition: monomials.h:44
EXTERN_VAR size_t gmp_output_digits
Definition: mpr_base.h:115
char * complexToStr(gmp_complex &c, const unsigned int oprec, const coeffs src)
Definition: mpr_complex.cc:704
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:60
#define nCopy(n)
Definition: numbers.h:15
#define nPrint(a)
only for debug, over any initalized currRing
Definition: numbers.h:46
#define nInit(i)
Definition: numbers.h:24
#define omFree(addr)
Definition: omAllocDecl.h:261
#define NULL
Definition: omList.c:12
static long pTotaldegree(poly p)
Definition: polys.h:282
#define pIsConstant(p)
like above, except that Comp must be 0
Definition: polys.h:238
#define pGetExp(p, i)
Exponent.
Definition: polys.h:41
void PrintS(const char *s)
Definition: reporter.cc:284
void PrintLn()
Definition: reporter.cc:310
static BOOLEAN rField_is_R(const ring r)
Definition: ring.h:519
static BOOLEAN rField_is_Zp(const ring r)
Definition: ring.h:501
static BOOLEAN rField_is_long_C(const ring r)
Definition: ring.h:546
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:507
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:593
@ LIST_CMD
Definition: tok.h:118
@ STRING_CMD
Definition: tok.h:185

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

4655{
4656 ideal gls = (ideal)(arg1->Data());
4657 int imtype= (int)(long)arg2->Data();
4658
4659 uResultant::resMatType mtype= determineMType( imtype );
4660
4661 // check input ideal ( = polynomial system )
4662 if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk )
4663 {
4664 return TRUE;
4665 }
4666
4667 uResultant *resMat= new uResultant( gls, mtype, false );
4668 if (resMat!=NULL)
4669 {
4670 res->rtyp = MODUL_CMD;
4671 res->data= (void*)resMat->accessResMat()->getMatrix();
4672 if (!errorreported) delete resMat;
4673 }
4674 return errorreported;
4675}
virtual ideal getMatrix()
Definition: mpr_base.h:31
const char * Name()
Definition: subexpr.h:120
Base class for solving 0-dim poly systems using u-resultant.
Definition: mpr_base.h:63
resMatrixBase * accessResMat()
Definition: mpr_base.h:78
VAR short errorreported
Definition: feFopen.cc:23
@ MODUL_CMD
Definition: grammar.cc:287
@ mprOk
Definition: mpr_base.h:98
uResultant::resMatType determineMType(int imtype)
mprState mprIdealCheck(const ideal theIdeal, const char *name, uResultant::resMatType mtype, BOOLEAN rmatrix=false)

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

4922{
4923 leftv v= args;
4924
4925 ideal gls;
4926 int imtype;
4927 int howclean;
4928
4929 // get ideal
4930 if ( v->Typ() != IDEAL_CMD )
4931 return TRUE;
4932 else gls= (ideal)(v->Data());
4933 v= v->next;
4934
4935 // get resultant matrix type to use (0,1)
4936 if ( v->Typ() != INT_CMD )
4937 return TRUE;
4938 else imtype= (int)(long)v->Data();
4939 v= v->next;
4940
4941 if (imtype==0)
4942 {
4943 ideal test_id=idInit(1,1);
4944 int j;
4945 for(j=IDELEMS(gls)-1;j>=0;j--)
4946 {
4947 if (gls->m[j]!=NULL)
4948 {
4949 test_id->m[0]=gls->m[j];
4950 intvec *dummy_w=id_QHomWeight(test_id, currRing);
4951 if (dummy_w!=NULL)
4952 {
4953 WerrorS("Newton polytope not of expected dimension");
4954 delete dummy_w;
4955 return TRUE;
4956 }
4957 }
4958 }
4959 }
4960
4961 // get and set precision in digits ( > 0 )
4962 if ( v->Typ() != INT_CMD )
4963 return TRUE;
4964 else if ( !(rField_is_R(currRing) || rField_is_long_R(currRing) || \
4966 {
4967 unsigned long int ii=(unsigned long int)v->Data();
4968 setGMPFloatDigits( ii, ii );
4969 }
4970 v= v->next;
4971
4972 // get interpolation steps (0,1,2)
4973 if ( v->Typ() != INT_CMD )
4974 return TRUE;
4975 else howclean= (int)(long)v->Data();
4976
4977 uResultant::resMatType mtype= determineMType( imtype );
4978 int i,count;
4979 lists listofroots= NULL;
4980 number smv= NULL;
4981 BOOLEAN interpolate_det= (mtype==uResultant::denseResMat)?TRUE:FALSE;
4982
4983 //emptylist= (lists)omAlloc( sizeof(slists) );
4984 //emptylist->Init( 0 );
4985
4986 //res->rtyp = LIST_CMD;
4987 //res->data= (void *)emptylist;
4988
4989 // check input ideal ( = polynomial system )
4990 if ( mprIdealCheck( gls, args->Name(), mtype ) != mprOk )
4991 {
4992 return TRUE;
4993 }
4994
4995 uResultant * ures;
4996 rootContainer ** iproots;
4997 rootContainer ** muiproots;
4998 rootArranger * arranger;
4999
5000 // main task 1: setup of resultant matrix
5001 ures= new uResultant( gls, mtype );
5002 if ( ures->accessResMat()->initState() != resMatrixBase::ready )
5003 {
5004 WerrorS("Error occurred during matrix setup!");
5005 return TRUE;
5006 }
5007
5008 // if dense resultant, check if minor nonsingular
5009 if ( mtype == uResultant::denseResMat )
5010 {
5011 smv= ures->accessResMat()->getSubDet();
5012#ifdef mprDEBUG_PROT
5013 PrintS("// Determinant of submatrix: ");nPrint(smv);PrintLn();
5014#endif
5015 if ( nIsZero(smv) )
5016 {
5017 WerrorS("Unsuitable input ideal: Minor of resultant matrix is singular!");
5018 return TRUE;
5019 }
5020 }
5021
5022 // main task 2: Interpolate specialized resultant polynomials
5023 if ( interpolate_det )
5024 iproots= ures->interpolateDenseSP( false, smv );
5025 else
5026 iproots= ures->specializeInU( false, smv );
5027
5028 // main task 3: Interpolate specialized resultant polynomials
5029 if ( interpolate_det )
5030 muiproots= ures->interpolateDenseSP( true, smv );
5031 else
5032 muiproots= ures->specializeInU( true, smv );
5033
5034#ifdef mprDEBUG_PROT
5035 int c= iproots[0]->getAnzElems();
5036 for (i=0; i < c; i++) pWrite(iproots[i]->getPoly());
5037 c= muiproots[0]->getAnzElems();
5038 for (i=0; i < c; i++) pWrite(muiproots[i]->getPoly());
5039#endif
5040
5041 // main task 4: Compute roots of specialized polys and match them up
5042 arranger= new rootArranger( iproots, muiproots, howclean );
5043 arranger->solve_all();
5044
5045 // get list of roots
5046 if ( arranger->success() )
5047 {
5048 arranger->arrange();
5049 listofroots= listOfRoots(arranger, gmp_output_digits );
5050 }
5051 else
5052 {
5053 WerrorS("Solver was unable to find any roots!");
5054 return TRUE;
5055 }
5056
5057 // free everything
5058 count= iproots[0]->getAnzElems();
5059 for (i=0; i < count; i++) delete iproots[i];
5060 omFreeSize( (ADDRESS) iproots, count * sizeof(rootContainer*) );
5061 count= muiproots[0]->getAnzElems();
5062 for (i=0; i < count; i++) delete muiproots[i];
5063 omFreeSize( (ADDRESS) muiproots, count * sizeof(rootContainer*) );
5064
5065 delete ures;
5066 delete arranger;
5067 if (smv!=NULL) nDelete( &smv );
5068
5069 res->data= (void *)listofroots;
5070
5071 //emptylist->Clean();
5072 // omFreeSize( (ADDRESS) emptylist, sizeof(slists) );
5073
5074 return FALSE;
5075}
int BOOLEAN
Definition: auxiliary.h:87
void * ADDRESS
Definition: auxiliary.h:119
Definition: intvec.h:23
virtual number getSubDet()
Definition: mpr_base.h:37
virtual IStateType initState() const
Definition: mpr_base.h:41
void solve_all()
Definition: mpr_numeric.cc:858
bool success()
Definition: mpr_numeric.h:162
void arrange()
Definition: mpr_numeric.cc:883
int getAnzElems()
Definition: mpr_numeric.h:95
rootContainer ** specializeInU(BOOLEAN matchUp=false, const number subDetVal=NULL)
Definition: mpr_base.cc:3059
rootContainer ** interpolateDenseSP(BOOLEAN matchUp=false, const number subDetVal=NULL)
Definition: mpr_base.cc:2921
@ denseResMat
Definition: mpr_base.h:65
@ IDEAL_CMD
Definition: grammar.cc:284
lists listOfRoots(rootArranger *self, const unsigned int oprec)
Definition: ipshell.cc:5078
#define nDelete(n)
Definition: numbers.h:16
#define nIsZero(n)
Definition: numbers.h:19
#define omFreeSize(addr, size)
Definition: omAllocDecl.h:260
void pWrite(poly p)
Definition: polys.h:308
int status int void size_t count
Definition: si_signals.h:59
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:35
intvec * id_QHomWeight(ideal id, const ring r)
#define IDELEMS(i)
Definition: simpleideals.h:23

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

4821{
4822 int i;
4823 ideal p,w;
4824 p= (ideal)arg1->Data();
4825 w= (ideal)arg2->Data();
4826
4827 // w[0] = f(p^0)
4828 // w[1] = f(p^1)
4829 // ...
4830 // p can be a vector of numbers (multivariate polynom)
4831 // or one number (univariate polynom)
4832 // tdg = deg(f)
4833
4834 int n= IDELEMS( p );
4835 int m= IDELEMS( w );
4836 int tdg= (int)(long)arg3->Data();
4837
4838 res->data= (void*)NULL;
4839
4840 // check the input
4841 if ( tdg < 1 )
4842 {
4843 WerrorS("Last input parameter must be > 0!");
4844 return TRUE;
4845 }
4846 if ( n != rVar(currRing) )
4847 {
4848 Werror("Size of first input ideal must be equal to %d!",rVar(currRing));
4849 return TRUE;
4850 }
4851 if ( m != (int)pow((double)tdg+1,(double)n) )
4852 {
4853 Werror("Size of second input ideal must be equal to %d!",
4854 (int)pow((double)tdg+1,(double)n));
4855 return TRUE;
4856 }
4857 if ( !(rField_is_Q(currRing) /* ||
4858 rField_is_R() || rField_is_long_R() ||
4859 rField_is_long_C()*/ ) )
4860 {
4861 WerrorS("Ground field not implemented!");
4862 return TRUE;
4863 }
4864
4865 number tmp;
4866 number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
4867 for ( i= 0; i < n; i++ )
4868 {
4869 pevpoint[i]=nInit(0);
4870 if ( (p->m)[i] )
4871 {
4872 tmp = pGetCoeff( (p->m)[i] );
4873 if ( nIsZero(tmp) || nIsOne(tmp) || nIsMOne(tmp) )
4874 {
4875 omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4876 WerrorS("Elements of first input ideal must not be equal to -1, 0, 1!");
4877 return TRUE;
4878 }
4879 } else tmp= NULL;
4880 if ( !nIsZero(tmp) )
4881 {
4882 if ( !pIsConstant((p->m)[i]))
4883 {
4884 omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4885 WerrorS("Elements of first input ideal must be numbers!");
4886 return TRUE;
4887 }
4888 pevpoint[i]= nCopy( tmp );
4889 }
4890 }
4891
4892 number *wresults= (number *)omAlloc( m * sizeof( number ) );
4893 for ( i= 0; i < m; i++ )
4894 {
4895 wresults[i]= nInit(0);
4896 if ( (w->m)[i] && !nIsZero(pGetCoeff((w->m)[i])) )
4897 {
4898 if ( !pIsConstant((w->m)[i]))
4899 {
4900 omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4901 omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4902 WerrorS("Elements of second input ideal must be numbers!");
4903 return TRUE;
4904 }
4905 wresults[i]= nCopy(pGetCoeff((w->m)[i]));
4906 }
4907 }
4908
4909 vandermonde vm( m, n, tdg, pevpoint, FALSE );
4910 number *ncpoly= vm.interpolateDense( wresults );
4911 // do not free ncpoly[]!!
4912 poly rpoly= vm.numvec2poly( ncpoly );
4913
4914 omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4915 omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4916
4917 res->data= (void*)rpoly;
4918 return FALSE;
4919}
Rational pow(const Rational &a, int e)
Definition: GMPrat.cc:411
int p
Definition: cfModGcd.cc:4078
vandermonde system solver for interpolating polynomials from their values
Definition: mpr_numeric.h:29
const CanonicalForm & w
Definition: facAbsFact.cc:51
#define nIsMOne(n)
Definition: numbers.h:26
#define nIsOne(n)
Definition: numbers.h:25
void Werror(const char *fmt,...)
Definition: reporter.cc:189