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

4648{
4649 res->data= (void*)loNewtonPolytope( (ideal)arg1->Data() );
4650 return FALSE;
4651}
#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 4653 of file ipshell.cc.

4654{
4655 if ( !(rField_is_long_R(currRing)) )
4656 {
4657 WerrorS("Ground field not implemented!");
4658 return TRUE;
4659 }
4660
4661 simplex * LP;
4662 matrix m;
4663
4664 leftv v= args;
4665 if ( v->Typ() != MATRIX_CMD ) // 1: matrix
4666 return TRUE;
4667 else
4668 m= (matrix)(v->CopyD());
4669
4670 LP = new simplex(MATROWS(m),MATCOLS(m));
4671 LP->mapFromMatrix(m);
4672
4673 v= v->next;
4674 if ( v->Typ() != INT_CMD ) // 2: m = number of constraints
4675 return TRUE;
4676 else
4677 LP->m= (int)(long)(v->Data());
4678
4679 v= v->next;
4680 if ( v->Typ() != INT_CMD ) // 3: n = number of variables
4681 return TRUE;
4682 else
4683 LP->n= (int)(long)(v->Data());
4684
4685 v= v->next;
4686 if ( v->Typ() != INT_CMD ) // 4: m1 = number of <= constraints
4687 return TRUE;
4688 else
4689 LP->m1= (int)(long)(v->Data());
4690
4691 v= v->next;
4692 if ( v->Typ() != INT_CMD ) // 5: m2 = number of >= constraints
4693 return TRUE;
4694 else
4695 LP->m2= (int)(long)(v->Data());
4696
4697 v= v->next;
4698 if ( v->Typ() != INT_CMD ) // 6: m3 = number of == constraints
4699 return TRUE;
4700 else
4701 LP->m3= (int)(long)(v->Data());
4702
4703#ifdef mprDEBUG_PROT
4704 Print("m (constraints) %d\n",LP->m);
4705 Print("n (columns) %d\n",LP->n);
4706 Print("m1 (<=) %d\n",LP->m1);
4707 Print("m2 (>=) %d\n",LP->m2);
4708 Print("m3 (==) %d\n",LP->m3);
4709#endif
4710
4711 LP->compute();
4712
4713 lists lres= (lists)omAlloc( sizeof(slists) );
4714 lres->Init( 6 );
4715
4716 lres->m[0].rtyp= MATRIX_CMD; // output matrix
4717 lres->m[0].data=(void*)LP->mapToMatrix(m);
4718
4719 lres->m[1].rtyp= INT_CMD; // found a solution?
4720 lres->m[1].data=(void*)(long)LP->icase;
4721
4722 lres->m[2].rtyp= INTVEC_CMD;
4723 lres->m[2].data=(void*)LP->posvToIV();
4724
4725 lres->m[3].rtyp= INTVEC_CMD;
4726 lres->m[3].data=(void*)LP->zrovToIV();
4727
4728 lres->m[4].rtyp= INT_CMD;
4729 lres->m[4].data=(void*)(long)LP->m;
4730
4731 lres->m[5].rtyp= INT_CMD;
4732 lres->m[5].data=(void*)(long)LP->n;
4733
4734 res->data= (void*)lres;
4735
4736 return FALSE;
4737}
#define TRUE
Definition: auxiliary.h:100
int m
Definition: cfEzgcd.cc:128
Variable next() const
Definition: factory.h:153
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:544
@ 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 4762 of file ipshell.cc.

4763{
4764 poly gls;
4765 gls= (poly)(arg1->Data());
4766 int howclean= (int)(long)arg3->Data();
4767
4768 if ( gls == NULL || pIsConstant( gls ) )
4769 {
4770 WerrorS("Input polynomial is constant!");
4771 return TRUE;
4772 }
4773
4775 {
4776 int* r=Zp_roots(gls, currRing);
4777 lists rlist;
4778 rlist= (lists)omAlloc( sizeof(slists) );
4779 rlist->Init( r[0] );
4780 for(int i=r[0];i>0;i--)
4781 {
4782 rlist->m[i-1].data=n_Init(r[i],currRing);
4783 rlist->m[i-1].rtyp=NUMBER_CMD;
4784 }
4785 omFree(r);
4786 res->data=rlist;
4787 res->rtyp= LIST_CMD;
4788 return FALSE;
4789 }
4790 if ( !(rField_is_R(currRing) ||
4794 {
4795 WerrorS("Ground field not implemented!");
4796 return TRUE;
4797 }
4798
4801 {
4802 unsigned long int ii = (unsigned long int)arg2->Data();
4803 setGMPFloatDigits( ii, ii );
4804 }
4805
4806 int ldummy;
4807 int deg= currRing->pLDeg( gls, &ldummy, currRing );
4808 int i,vpos=0;
4809 poly piter;
4810 lists elist;
4811
4812 elist= (lists)omAlloc( sizeof(slists) );
4813 elist->Init( 0 );
4814
4815 if ( rVar(currRing) > 1 )
4816 {
4817 piter= gls;
4818 for ( i= 1; i <= rVar(currRing); i++ )
4819 if ( pGetExp( piter, i ) )
4820 {
4821 vpos= i;
4822 break;
4823 }
4824 while ( piter )
4825 {
4826 for ( i= 1; i <= rVar(currRing); i++ )
4827 if ( (vpos != i) && (pGetExp( piter, i ) != 0) )
4828 {
4829 WerrorS("The input polynomial must be univariate!");
4830 return TRUE;
4831 }
4832 pIter( piter );
4833 }
4834 }
4835
4836 rootContainer * roots= new rootContainer();
4837 number * pcoeffs= (number *)omAlloc( (deg+1) * sizeof( number ) );
4838 piter= gls;
4839 for ( i= deg; i >= 0; i-- )
4840 {
4841 if ( piter && pTotaldegree(piter) == i )
4842 {
4843 pcoeffs[i]= nCopy( pGetCoeff( piter ) );
4844 //nPrint( pcoeffs[i] );PrintS(" ");
4845 pIter( piter );
4846 }
4847 else
4848 {
4849 pcoeffs[i]= nInit(0);
4850 }
4851 }
4852
4853#ifdef mprDEBUG_PROT
4854 for (i=deg; i >= 0; i--)
4855 {
4856 nPrint( pcoeffs[i] );PrintS(" ");
4857 }
4858 PrintLn();
4859#endif
4860
4861 roots->fillContainer( pcoeffs, NULL, 1, deg, rootContainer::onepoly, 1 );
4862 roots->solver( howclean );
4863
4864 int elem= roots->getAnzRoots();
4865 char *dummy;
4866 int j;
4867
4868 lists rlist;
4869 rlist= (lists)omAlloc( sizeof(slists) );
4870 rlist->Init( elem );
4871
4873 {
4874 for ( j= 0; j < elem; j++ )
4875 {
4876 rlist->m[j].rtyp=NUMBER_CMD;
4877 rlist->m[j].data=(void *)nCopy((number)(roots->getRoot(j)));
4878 //rlist->m[j].data=(void *)(number)(roots->getRoot(j));
4879 }
4880 }
4881 else
4882 {
4883 for ( j= 0; j < elem; j++ )
4884 {
4885 dummy = complexToStr( (*roots)[j], gmp_output_digits, currRing->cf );
4886 rlist->m[j].rtyp=STRING_CMD;
4887 rlist->m[j].data=(void *)dummy;
4888 }
4889 }
4890
4891 elist->Clean();
4892 //omFreeSize( (ADDRESS) elist, sizeof(slists) );
4893
4894 // this is (via fillContainer) the same data as in root
4895 //for ( i= deg; i >= 0; i-- ) nDelete( &pcoeffs[i] );
4896 //omFreeSize( (ADDRESS) pcoeffs, (deg+1) * sizeof( number ) );
4897
4898 delete roots;
4899
4900 res->data= (void*)rlist;
4901
4902 return FALSE;
4903}
#define NULL
Definition: auxiliary.h:104
int i
Definition: cfEzgcd.cc:132
int * Zp_roots(const CanonicalForm f)
Definition: cf_roots.cc:25
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:299
int getAnzRoots()
Definition: mpr_numeric.h:97
bool solver(const int polishmode=PM_NONE)
Definition: mpr_numeric.cc:436
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:539
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
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:520
static BOOLEAN rField_is_Zp(const ring r)
Definition: ring.h:502
static BOOLEAN rField_is_long_C(const ring r)
Definition: ring.h:547
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:508
static short rVar(const ring r)
#define rVar(r) (r->N)
Definition: ring.h:594
@ 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 4739 of file ipshell.cc.

4740{
4741 ideal gls = (ideal)(arg1->Data());
4742 int imtype= (int)(long)arg2->Data();
4743
4744 uResultant::resMatType mtype= determineMType( imtype );
4745
4746 // check input ideal ( = polynomial system )
4747 if ( mprIdealCheck( gls, arg1->Name(), mtype, true ) != mprOk )
4748 {
4749 return TRUE;
4750 }
4751
4752 uResultant *resMat= new uResultant( gls, mtype, false );
4753 if (resMat!=NULL)
4754 {
4755 res->rtyp = MODUL_CMD;
4756 res->data= (void*)resMat->accessResMat()->getMatrix();
4757 if (!errorreported) delete resMat;
4758 }
4759 return errorreported;
4760}
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 5006 of file ipshell.cc.

5007{
5008 leftv v= args;
5009
5010 ideal gls;
5011 int imtype;
5012 int howclean;
5013
5014 // get ideal
5015 if ( v->Typ() != IDEAL_CMD )
5016 return TRUE;
5017 else gls= (ideal)(v->Data());
5018 v= v->next;
5019
5020 // get resultant matrix type to use (0,1)
5021 if ( v->Typ() != INT_CMD )
5022 return TRUE;
5023 else imtype= (int)(long)v->Data();
5024 v= v->next;
5025
5026 if (imtype==0)
5027 {
5028 ideal test_id=idInit(1,1);
5029 int j;
5030 for(j=IDELEMS(gls)-1;j>=0;j--)
5031 {
5032 if (gls->m[j]!=NULL)
5033 {
5034 test_id->m[0]=gls->m[j];
5035 intvec *dummy_w=id_QHomWeight(test_id, currRing);
5036 if (dummy_w!=NULL)
5037 {
5038 WerrorS("Newton polytope not of expected dimension");
5039 delete dummy_w;
5040 return TRUE;
5041 }
5042 }
5043 }
5044 }
5045
5046 // get and set precision in digits ( > 0 )
5047 if ( v->Typ() != INT_CMD )
5048 return TRUE;
5049 else if ( !(rField_is_R(currRing) || rField_is_long_R(currRing) || \
5051 {
5052 unsigned long int ii=(unsigned long int)v->Data();
5053 setGMPFloatDigits( ii, ii );
5054 }
5055 v= v->next;
5056
5057 // get interpolation steps (0,1,2)
5058 if ( v->Typ() != INT_CMD )
5059 return TRUE;
5060 else howclean= (int)(long)v->Data();
5061
5062 uResultant::resMatType mtype= determineMType( imtype );
5063 int i,count;
5064 lists listofroots= NULL;
5065 number smv= NULL;
5066 BOOLEAN interpolate_det= (mtype==uResultant::denseResMat)?TRUE:FALSE;
5067
5068 //emptylist= (lists)omAlloc( sizeof(slists) );
5069 //emptylist->Init( 0 );
5070
5071 //res->rtyp = LIST_CMD;
5072 //res->data= (void *)emptylist;
5073
5074 // check input ideal ( = polynomial system )
5075 if ( mprIdealCheck( gls, args->Name(), mtype ) != mprOk )
5076 {
5077 return TRUE;
5078 }
5079
5080 uResultant * ures;
5081 rootContainer ** iproots;
5082 rootContainer ** muiproots;
5083 rootArranger * arranger;
5084
5085 // main task 1: setup of resultant matrix
5086 ures= new uResultant( gls, mtype );
5087 if ( ures->accessResMat()->initState() != resMatrixBase::ready )
5088 {
5089 WerrorS("Error occurred during matrix setup!");
5090 return TRUE;
5091 }
5092
5093 // if dense resultant, check if minor nonsingular
5094 if ( mtype == uResultant::denseResMat )
5095 {
5096 smv= ures->accessResMat()->getSubDet();
5097#ifdef mprDEBUG_PROT
5098 PrintS("// Determinant of submatrix: ");nPrint(smv);PrintLn();
5099#endif
5100 if ( nIsZero(smv) )
5101 {
5102 WerrorS("Unsuitable input ideal: Minor of resultant matrix is singular!");
5103 return TRUE;
5104 }
5105 }
5106
5107 // main task 2: Interpolate specialized resultant polynomials
5108 if ( interpolate_det )
5109 iproots= ures->interpolateDenseSP( false, smv );
5110 else
5111 iproots= ures->specializeInU( false, smv );
5112
5113 // main task 3: Interpolate specialized resultant polynomials
5114 if ( interpolate_det )
5115 muiproots= ures->interpolateDenseSP( true, smv );
5116 else
5117 muiproots= ures->specializeInU( true, smv );
5118
5119#ifdef mprDEBUG_PROT
5120 int c= iproots[0]->getAnzElems();
5121 for (i=0; i < c; i++) pWrite(iproots[i]->getPoly());
5122 c= muiproots[0]->getAnzElems();
5123 for (i=0; i < c; i++) pWrite(muiproots[i]->getPoly());
5124#endif
5125
5126 // main task 4: Compute roots of specialized polys and match them up
5127 arranger= new rootArranger( iproots, muiproots, howclean );
5128 arranger->solve_all();
5129
5130 // get list of roots
5131 if ( arranger->success() )
5132 {
5133 arranger->arrange();
5134 listofroots= listOfRoots(arranger, gmp_output_digits );
5135 }
5136 else
5137 {
5138 WerrorS("Solver was unable to find any roots!");
5139 return TRUE;
5140 }
5141
5142 // free everything
5143 count= iproots[0]->getAnzElems();
5144 for (i=0; i < count; i++) delete iproots[i];
5145 omFreeSize( (ADDRESS) iproots, count * sizeof(rootContainer*) );
5146 count= muiproots[0]->getAnzElems();
5147 for (i=0; i < count; i++) delete muiproots[i];
5148 omFreeSize( (ADDRESS) muiproots, count * sizeof(rootContainer*) );
5149
5150 delete ures;
5151 delete arranger;
5152 nDelete( &smv );
5153
5154 res->data= (void *)listofroots;
5155
5156 //emptylist->Clean();
5157 // omFreeSize( (ADDRESS) emptylist, sizeof(slists) );
5158
5159 return FALSE;
5160}
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:857
bool success()
Definition: mpr_numeric.h:162
void arrange()
Definition: mpr_numeric.cc:882
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:5163
#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 4905 of file ipshell.cc.

4906{
4907 int i;
4908 ideal p,w;
4909 p= (ideal)arg1->Data();
4910 w= (ideal)arg2->Data();
4911
4912 // w[0] = f(p^0)
4913 // w[1] = f(p^1)
4914 // ...
4915 // p can be a vector of numbers (multivariate polynom)
4916 // or one number (univariate polynom)
4917 // tdg = deg(f)
4918
4919 int n= IDELEMS( p );
4920 int m= IDELEMS( w );
4921 int tdg= (int)(long)arg3->Data();
4922
4923 res->data= (void*)NULL;
4924
4925 // check the input
4926 if ( tdg < 1 )
4927 {
4928 WerrorS("Last input parameter must be > 0!");
4929 return TRUE;
4930 }
4931 if ( n != rVar(currRing) )
4932 {
4933 Werror("Size of first input ideal must be equal to %d!",rVar(currRing));
4934 return TRUE;
4935 }
4936 if ( m != (int)pow((double)tdg+1,(double)n) )
4937 {
4938 Werror("Size of second input ideal must be equal to %d!",
4939 (int)pow((double)tdg+1,(double)n));
4940 return TRUE;
4941 }
4942 if ( !(rField_is_Q(currRing) /* ||
4943 rField_is_R() || rField_is_long_R() ||
4944 rField_is_long_C()*/ ) )
4945 {
4946 WerrorS("Ground field not implemented!");
4947 return TRUE;
4948 }
4949
4950 number tmp;
4951 number *pevpoint= (number *)omAlloc( n * sizeof( number ) );
4952 for ( i= 0; i < n; i++ )
4953 {
4954 pevpoint[i]=nInit(0);
4955 if ( (p->m)[i] )
4956 {
4957 tmp = pGetCoeff( (p->m)[i] );
4958 if ( nIsZero(tmp) || nIsOne(tmp) || nIsMOne(tmp) )
4959 {
4960 omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4961 WerrorS("Elements of first input ideal must not be equal to -1, 0, 1!");
4962 return TRUE;
4963 }
4964 } else tmp= NULL;
4965 if ( !nIsZero(tmp) )
4966 {
4967 if ( !pIsConstant((p->m)[i]))
4968 {
4969 omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4970 WerrorS("Elements of first input ideal must be numbers!");
4971 return TRUE;
4972 }
4973 pevpoint[i]= nCopy( tmp );
4974 }
4975 }
4976
4977 number *wresults= (number *)omAlloc( m * sizeof( number ) );
4978 for ( i= 0; i < m; i++ )
4979 {
4980 wresults[i]= nInit(0);
4981 if ( (w->m)[i] && !nIsZero(pGetCoeff((w->m)[i])) )
4982 {
4983 if ( !pIsConstant((w->m)[i]))
4984 {
4985 omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
4986 omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
4987 WerrorS("Elements of second input ideal must be numbers!");
4988 return TRUE;
4989 }
4990 wresults[i]= nCopy(pGetCoeff((w->m)[i]));
4991 }
4992 }
4993
4994 vandermonde vm( m, n, tdg, pevpoint, FALSE );
4995 number *ncpoly= vm.interpolateDense( wresults );
4996 // do not free ncpoly[]!!
4997 poly rpoly= vm.numvec2poly( ncpoly );
4998
4999 omFreeSize( (ADDRESS)pevpoint, n * sizeof( number ) );
5000 omFreeSize( (ADDRESS)wresults, m * sizeof( number ) );
5001
5002 res->data= (void*)rpoly;
5003 return FALSE;
5004}
Rational pow(const Rational &a, int e)
Definition: GMPrat.cc:411
int p
Definition: cfModGcd.cc:4080
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